Theory of double-resonant Raman spectra in graphene: intensity and line shape of 

defect-induced and two-phonon bands 
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We calculate the double resonant (DR) Raman spectrum of graphene, and determine the lines 
associated to both phonon-defect processes (such as in the D line at ~ 1350 cm" 1 , D' at ~ 1600 cm" 
and D" at ~ 1100 cm -1 ), and two-phonons ones (such as in the 2D, 2D', or D + D" lines). Phonon 
and electronic dispersions reproduce calculations based on density functional theory corrected with 
GW. Electron-light, -phonon , and -defect scattering matrix elements and the electronic linewidth 
are explicitly calculated. Defect-induced processes are simulated by considering different kind of 
idealized defects. For an excitation energy of el = 2.4 eV, the agreement with measurements is very 
good and calculations reproduce: the relative intensities among phonon-defect or among two-phonon 
lines; the measured small widths of the D, D' , 2D and 2D' lines; the line shapes; the presence of 
small intensity lines in the 1800, 2000 cm -1 range. We determine how the spectra depend on the 
excitation energy, on the light polarization, on the electronic linewidth, on the kind of defects and 
on their concentration. According to the present findings, the intensity ratio between the 2D' and 
2D lines can be used to determine experimentally the electronic linewidth. The intensity ratio 
between the D and D' lines depends on the kind of model defect, suggesting that this ratio could 
possibly be used to identify the kind of defects present in actual samples. Charged impurities 
outside the graphene plane provide an almost undetectable contribution to the Raman signal. The 
present analysis reveals that, for both D and 2D lines, the dominant DR processes are those in 
which electrons and holes are both involved in the scattering, because of a destructive quantum 
interference that kills processes involving only electrons or only holes. The most important phonons 
belong to the K— s> T direction (inner phonons) and not to the K— >M one (outer phonons), as 
usually assumed. The small 2D line width at cl — 2.4 eV is a consequence of the interplay between 
the opposite trigonal warpings of the electron and phonon dispersions. At higher excitation, e.g. 
e_L = 3.8 eV, the 2D line becomes broader and evolves in an asymmetric double peak structure. 

PACS numbers: 78.30.-j,78.67.Wj,81.05.ue 



I. INTRODUCTION 

Raman spectroscopy is one of the most impor- 
tant experimental techniques for the characterization of 
graphitic materials. In particular, for graphene, this 
technique provides information about the number of lay- 
ers [11,0) doping [HH, disorder and phonon prop- 
erties Q. 

Lowest-oder Raman processes correspond to the scat- 
tering with a zero momentum phonon (q=0). The Ra- 
man G line in graphene and graphite (^1582 cm" 1 ) is as- 
sociated with the E2 9 phonon at T and it is a lowest-order 
process. Graphene and graphite present other lines, due 
to higher order processes, which are usually interpreted 
in terms of the so called double resonance (DR) mecha- 
nism [10]. The DR mechanism is used to interpret two 
distinct kind of phenomena. The first is the excitation 
of a phonon with momentum q^0 due to the presence of 
defects in the sample. This process, called defect-induced, 
is not allowed in a purely crystalline sample (without de- 
fects) because of momentum conservation. In graphene 
and graphite, it gives rise to the well studied D line at 
~ 1350 cm" 1 and also to less intense lines such as the 
D' (~ 1600 cm" 1 ), and the D" (~ 1100 cm" 1 @, [lH ) 
The second process corresponds to the excitation of two 
phonons with opposite momenta q and -q. This process, 
called two-phonon, can be observed in purely crystalline 



samples since the momentum is conserved and gives rise 
to the very intense 2D line at ~ 2700 cm" 1 (which is an 
overtone of the D line) and, for instance, to the D + D" 
and 2D' lines at - 2450 cm" 1 and - 3200 cm" 1 . The 
lines related to DR defect-induced and two-phonon pro- 
cesses have a remarkable property: they are dispersive, 
i.e. their positions change with excitation energy. 

It has been shown experimentally [j], [|[ that the 2D 
line in graphene changes in shape, width and position 
with number of layers. Later, the phonon dispersion of 
graphene, near the Dirac K points, was probed by mea- 
surements @ of the 2D and D + D" lines as a function 
of the excitation energies. Usually, Raman experiments 
are performed in graphene layers that were deposited or 
grown over a substrate. However, experimental measure- 
ments of the G and 2D lines have also been performed 
for free-standing graphene monolayers 12j . Lucchese et. 
al 0] and Martins Ferreira et. al [ll[ have studied the 
evolution of the Raman spectra for mono and multi-layer 
graphene with increasing disorder, showing that the in- 
tensity of the D line, which is absent in pristine graphene, 
increases when disorder is induced in the sample up to 
a maximum value where it begins to decrease. On the 
other hand, the 2D line intensity is maximum for pristine 
graphene and it decreases with increasing disorder. 

Frequencies, intensities and linewidths of all DR Ra- 
man bands may be determined by the calculation of the 
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Raman cross section Several excellent theoretical 

works already appeared on the topic providing an over- 
all good understanding of the situation. However, the 
many different approximations used by different authors 
(e.g. constant electron-phonon matrix elements, reso- 
nant phonons are assumed to be on some high symmetry 
line, in some cases the electronic dispersion is conic, the 
electronic life-time is a parameter, etc.) and the sev- 
eral debates still going on lead the sensation that some- 
thing is missing. Thomsen and Reich fiolj and Kurti et. 
al [lH studied the D line for graphite and carbon nan- 
otubes, respectively. Also, Narula and Reich[lj| studied 
the D and 2D Raman lines in graphene and graphite. In 
these works pH [HI, HH the scattering matrix elements 
(electron-light, electron-phonon and electron-defect) are 
assumed to be constants and the electronic linewidth is 
a parameter set ot a fixed value. Basko [l6[ has stud- 
ied the two-phonon and four-phonon Raman bands in 
graphene under the assumption of conical bands, which 
is valid only in the limit of small excitation energies, not 
suitable for most experimental data available in the liter- 
ature. Also, his work is limited to disorder-free graphene. 
Park et. al [l7| have studied the two-phonon processes 
in single, double and triple layer graphene, making the 
assumption of conical bands and limiting their work to 
disorder-free graphene. 

In this context, some questions are currently de- 
bated. For instance, according to previous theoretical 
works [l(J EH, Ell j phonons in the K— >M direction of the 
Brillouin zone should give the most important contribu- 
tion to the D line intensity. However, recent works p~8l — 
I22I ] have argued that the phonons in the K— ^ T direction 
should be more important. Other open questions refer to 
the processes more relevant for the DR Raman spectra. 
In some Raman processes only the electrons are scat- 
tered, while in other processes both electrons and holes 
are scattered simultaneously. Some authors claim that, 
at least for the 2D line, this last kind of processes should 
be dominant because they are associated to a triple res- 
onance (23[. On the other hand, several authors perform 



their studies considering only electron-electron processes, 
as in the seminal work by Thomsen and Reich ■ 

Besides, several fundamental questions are almost un- 
touched. So far, the DR mechanism has been basically 
used to give an overall description of the physics and to 
determine which are the excited phonons. Can the DR 
theory be used to obtain a quantitative description of 
the intensities of the Raman lines? Can the DR theory 
be used to obtain a quantitative description of the shape 
and of the width of the Raman lines? The most studied 
Raman lines, the 2D and the D ones, present a relatively 
narrow linewidth similar to the one of the G line (which is 
not due to DR). This fact is very surprising and, indeed, 
the theoretical approaches used so far were not able to 
reproduce the observed small width of these lines. Which 
are the missing ingredients? Is this a consequence of the 
approximations used so far, or, on the contrary, is this 
a limit of the perturbative approach inherent to the DR 



theory? Finally, the D line is activated by disorder and 
is routinely used to probe the quality of the samples of 
graphitic materials. However, which kind of defects ac- 
tivate the D line is not known. For instance, do neutral 
impurities, vacancies and charged defects affect the D 
line in the same way? Which kind of defects are probed 
by measuring different defect- activates lines? Does Ra- 
man spectroscopy probe the defects which mostly influ- 
ence electronic transport? 

Here, as a first step to answer these questions, we cal- 
culate the double resonant Raman spectrum of graphene, 
considering both defect-induced and two-phonons pro- 
cesses, trying to provide a computational method over- 
coming the most common approximations used in litera- 
ture. Calculations are done using the standard approach 
based on the golden rule generalized to the perturba- 
tive fourth-order (To| . The electronic summation is per- 
formed all over the two dimensional Brillouin zone and 
all the possible phonons (with any wavevector) are con- 
sidered. The phonon dispersion is obtained from fully 
ab-initio calculations based on density functional theory 
(DFT) corrected with GW. Electronic structure calcula- 
tions are based on a tight binding approach in which the 
parameters are fitted to reproduce DFT+GW calcula- 
tions. The electronic lifetime is calculated explicitly and 
the defect-induced processes are simulated by considering 
three different kind of ideal model defects. 

Sec. |TT] describes the computational method; Sec. IIIII 
describes and discusses the results; Sec. IIVI resumes the 
main conclusions of the paper. 



II. METHOD 

This section describes the method used to compute the 
DR Raman spectra. Sec. Ill Al gives the general frame- 
work and provides the equations to obtain double res- 
onant Raman spectra in graphene within the perturba- 
tive approach. The other subsections describe the details 
to obtain the quantities used in the actual implemen- 
tation. In particular, Sec. Ill Bl describes the electronic 
and phononic band dispersions; Sects. HTCl III Dl iHEl dc- 
scribe the electron-phonon, electron-light and electron- 
defect scattering matrix elements; Sec. Ill Fl describes the 
calculation of the electronic linewidth. 



A. Double resonant Raman intensity 

In vibrational Raman, the spectrum usually consists 
in well defined lines associated with emission (Stokes) 
or absorption (anti-Stokcs) of a phonon. Here, only 
Stokes processes are considered. Note also that the G line 
(lowest-order excitation of the E2 ff T phonon) is not de- 
scribed by the present formalism and is, thus, not present 
in the calculated spectra. Within the DR scheme [Toj . 
the light-electron and electron-phonon interactions, as 
well as the defect-induced electron-electron scattering are 
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treated at the first order in perturbation theory. The Ra- obtained from the golden rule generalized to the fourth- 
man cross section / of the light scattered by a crystal is order [13j : 
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where is the energy of the initial state which consists 
in a quantum of light with energy — Hujl (the laser en- 
ergy) and in which the crystal is in the ground state. The 
sum is performed on intermediate virtual states A, B, C, 
with energy ca, c-b, ecs which are described by electronic 
and phononic excitation of the crystal, et is the energy 
of the final state /, in which the electronic degrees of 
freedom of the crystal are in the ground state, one or 
two phonons with total energy hu) p have been excited, 
and a quantum of light with energy ex, — Hu) p has been 
emitted. 8 is the Dirac distribution. 7 A , 7 s , j c are the 
inverse of the lifetimes of the electronic excitations of the 
virtual states A, B, C, respectively. M.jk are first-order 
scattering matrix elements between the states J and K. 
So far, no attempts have been reported to go beyond the 
approximation inherent to Eq. [U for graphitic materials. 
Note that within the present approach, the G line (which 
in literature is usually referred to as a "first-order" pro- 
cess) is a third-order process. 

The processes described by Eq. [T] are in general asso- 
ciated to lines which are much weaker than "first-order" 
Raman lines. Graphene and graphite are notable excep- 
tions. During the intermediate virtual transition the en- 
ergy is not necessarily conserved and the three denomina- 
tors of Eq. Q]are in general different from zero. However, 



in graphene and graphite two or more of the denomina- 
tors of Eq. [T] can be equal to zero simultaneously. In 
literature this is called double-resonance condition, and 
can be associated to Raman lines which have an intensity 
comparable to that of lower-order processes (the G line) . 

In the DR Raman scattering, the process Mai m Eq.Q] 
corresponds to the absorption of light by creation of an 
electron- hole pair in the tt/tt* bands. Then, the carri- 
ers are scattered twice before recombination {M ba and 
Mcb in Eq. [T]). For temperatures typically present in 
Raman measurements in graphene, only Stoke processes 
(phonon emission) are relevant. Thus, in one possible 
case, one scattering event is due to collision with a de- 
fect and the other to the creation of a phonon {phonon- 
defect process). In a second possible case, both scattering 
events are due to creation of phonons (two-phonon pro- 
cess). Finally, the process Ai/c in Eq. Q] corresponds to 
the recombination of the carriers by light emission. We 
define 1?% as the probability to excite a phonon -qy, with 
momentum -q, branch index v and energy Sw^. q through 
a phonon- defect process. Ifffia is the probability to excite 
the two phonons -qy and q/x through a two-phonon pro- 
cess. The Raman intensity as a function of the frequency 
u) of the scattered light is proportional to 
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The sum in Eq. [5] is performed on a uniform grid of N q 
phonon wavevectors q in the Brillouin zone and on all 
the branch indexes v and /i. In the limit iV q — > 00, S(uj) 
is the Dirac distribution. n(w) is the Bose-Einstein oc- 
cupation. In Eq. [21 Nd is the average number of defects 
in the unit cell. I pd cx Nd, because we assume that the 
contributions of defects on different sites add up incoher- 
ently. The first sum in Eq. [3] is performed on a uniform 



grid of Nk electronic wavevectors k. a and /3 are labels 
running on the eight different possible processes that we 
call eel, ee2, hhl, hh2, ehl, eh2, hel, he2, which are rep- 
resented diagrammatically in Fig. [TJ The reader might 
be familiar with an alternative representation of the pro- 
cesses, reported in Fig. [2] Expressions for the DR scat- 
tering amplitudes K are given in the appendix. Here we 
report, as examples, K^ el and K^: 
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Phonon-defect (pd) processes: 



eel 



-qv 



ee2 



hhl 



^-qv 



hh2 



-qv 





„-qv 


ehl 




eh2 






Sj <-qv 


hel 


» 




-Q 

s ^-qv 






he2 








eel 



ee2 



hhl 



hh2 



Light: 



Electron: 



Hole: 



Two-phonons (pp) processes: 

^-qv^qji 



? qn ^ -qv 



Jk-qv^q|j. 



* q|i A -qv 



Phonon: 





„-qv 




s *qH 


hel 


*-qv 




he2 


s ^-qv 

„-qv 


s *qn 


Defect: 


X 



FIG. 1: Goldstone diagrams for the double resonant Raman processes considered in this work. In this manuscript, the term 
"ai> processes" refers to the processes highlighted by the gray area (ehl, eh2, hel, and he2). The other processes are referred 
to as "aa processes". The largest part of the Raman intensity is due to the ab processes. The reader might be familiar with an 
alternative representation of the processes, reported in Fig. [21 
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Eq. [4] corresponds to the phonon-defect diagram eel in 
Fig. [TJ Initially, the excitation laser creates an electron- 
hole pair with momentum k. Thus, using the notation 
of Eq. [TJ Mai = (7r* k| Di n 1 7rk), where |k7r) and |k7r*) 
are the electronic occupied and empty states and is 
the operator coupling the incident electromagnetic wave 
with the crystal, ej = ez, and = el — si, being e£ 
the energy of |k7r). Secondly, the excited electron is scat- 
tered into a k + q state by emitting a phonon with mo- 
mentum -q. Thus, Mba = (k + q, 7r*|AiJ qj „|k7r*), be- 
ing Aff q l/ the electron-phonon coupling operator. Now, 



£b 



-el+hw^. 



The third step in the process K\ 



■pd 



is the scattering of the k + q electron by a defect back to 
the k state. Thus, Mcb = (k7r*|ifu|k+q, tt*), being H D 
the defect scattering operator and ec = el — el + hui^. 
Finally, the electron and hole recombine vertically in the 
k-state, by emitting light. Thus, Mfc — (k7r|£J out |k7r*), 
being D out the operator coupling the emitted photon 
with the crystal. The broadening energies 7k in the de- 
nominators of the DR amplitudes K (e.g. in Eqs.|4l [5]) 



are the inverse of the corresponding electronic lifetimes 
(see SecHm. 

Eq. O corresponds to the phonon-phonon diagram 
eel in Fig. [T] The first two step are the same as in 
the previous paragraph, while in the third step, the 
k + q electron is scattered into a k electron, by emit- 
ting the phonon with momentum q/^. Thus, A4cb = 
(kTr-IAiL-qJk + qir*) and e c = e£* - e£ + ^ q + 
The fourth step is the same as before. Finally, for 
graphene and graphite, the diagrams of Fig. [T]are some- 
times schematized with a different notation. For a com- 
parison see Fig. [21 

The sums in Eq. [2] are performed on a uniform grid 
of 120x120 q points (randomly shifted with respect to 
the origin) and 5{oS) is a Lorentzian distribution with 
8 cm -1 full width at half maximum. The results will be 
plotted as a function of the Raman shift ujl — The 
sums in Eq. [3l are performed on grids of k points which 
are sufficiently large to ensure convergence. Depending 
on the value of 7^ uniform grids between 480x480 and 
840x840 k points are used. In Eq.[3J we consider fiw q ^> 
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FIG. 2: An alternative representation (customary for 
graphene and graphite) of the processes associated to the di- 
agrams of Fig. [T] The crosses represent the electronic disper- 
sion near the conic region. The vertical arrows represent the 
electron/hole creation and recombination. The horizontal ar- 
rows represent the scattering with a defect or with a phonon. 
For simplicity we show only the processes involving a phonon 
with momentum along the K-M line. In this manuscript, the 
term u ab processes" refers to the processes highlighted by the 
gray area (ehl, eh2, hel, and he2). The other processes are 
referred to as "aa processes". 



tic which is badly described by small-neighbors TBs, but 
which is well described by the present 5-neighbors TB, 
is the electron/hole asymmetry, e£* + e£. This quantity 
depends on the k direction and has values of the order 
of the electronic broadening (see Sec. IIIFj) : e.g. for the 
states in resonance with a laser of 2.4 eV, the asymme- 
try is about 40 and 100 meV along the K~r and the 
K-M direction, respectively (Fig. [3]). On the contrary, 
in a lst-neighbors TB model, the e/h asymmetry is k 
independent and it is equal to zero. 



3 



-3 











vs. 

— 

~ // 
s / 




>^ — - nearest-neighbours 
^ — up to 5th neighbours - 





k m 



KbT and, thus, n(wq) ~ 0. Unless otherwise specified, 
the intensities are normalized to the maximum value of 
the 2D peak. In the following four sub-sections (and 
in App. [Bjl . we describe the model to obtain the DR 
scattering amplitudes K. 

B. Electron and phonon dispersion 

The electronic structure, e£ and |k, a), is obtained 
from a tight binding (TB) model with one orthonormal- 
ized p z orbital per site and interactions up to fifth neigh- 
bors (details are in App. IB II) . We use t\ = —3.40 eV, 
t 2 = 0.33 eV, t 3 = -0.24 eV, U = 0.12 eV and t 5 = 0.09 
eV, where U is the i-th neighbor hopping parameter. The 
resulting electronic dispersion is shown in Fig. [3J These 
TB parameters were obtained following (24|: first, the U 
are fitted to density-functional theory (DFT) electronic 
band dispersion to reproduce the tt — tt* bands along the 
T-K-M line; then, all the U are rescaled by +18% in 
order to reproduce the 7r band slope near K from GW 
calculations, which are in excellent agreement with angle- 
resolved photoemission spectra (ARPES) measurements 
on graphite [25j |. 

We remark that, in the present context, a good descrip- 
tion of the trigonal warping of the 7r-bands cone is very 
relevant, since the actual shape of the trigonal warping 
determines the q vectors of the phonons associated to the 
D line. The present 5-neighbors TB can reproduce very 
well the trigonal warping as obtained from DFT. On the 
contrary, by using a lst-neighbors TB model, the trigonal 
warping is underestimated. Another relevant characteris- 



FIG. 3: Graphene electronic dispersion obtained with the 
5-neighbors tight-binding described in the text (solid line). 
For a comparison we also show the dispersion obtained with 
the 1st neighbors TB having the same Fermi velocity at K 
(dashed line). The electron/hole (e/h) asymmetry is defined 
as e£ + e£ and is constant for the 1st neighbors TB. 

Phonon dispersions, w^, are obtained from ab-initio 
DFT calculations [26| corrected with GW as in [III, SI- 
In particular, first we computed the DFT phonon disper- 
sion, then we "correct" the dispersion of the highest op- 
tical branch near K (the branch which is TO near T and 
which is associated with the A^ mode at K, see Fig. 2]) 
by rescaling the phonon self-energy contribution to the 
dynamical matrix consistently with the GW calculated 
clcctron-phonon coupling and electronic ir band disper- 
sion [27]]. Calculations are done for graphene with the 
same computational details of [28( . In [28| , the rescaling 
factor is a constant, r GW = 1.61, all over the BZ and the 
phonons are studied just in the neighborhood of K. Here, 
in order to obtain a phonon dispersion all over the BZ, 
the rescaling factor, r GW , depends on q. r GW — r GW 
near K and smoothly drops to one elsewhere: 

r GW = 1 + (r GW - l)-erfc ( |q ~ " — ) , (6) 

q 1 '2 J \ 0.05 / 

being ao the graphene lattice constant and K n the near- 
est vector to q among those equivalent to K. The GW 
correction associated to r GW changes the phonon slope of 
the highest optical branch near K by almost +60% (with 
respect to DFT) providing a much better agreement with 
measurements for graphite (Fig. 0|) . The precise value of 
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the phonon dispersion near K is essential in the present 
context, since it determines the dependence of the D peak 
dispersion as a function of the exciting laser energy [2!| . 

Finally, notice that the present DFT calculations re- 
produce very well the experimental phonon dispersion 
from inelastic x-ray scattering (IXS) of [30| of the high- 
est optical branch near T. We can thus assume that the 
DFT frequency for the E2 9 T mode (1561 cm -1 ) is a pre- 
cise fit of the IXS measurements. The 1561 cm -1 value is 
however 1.3% smaller than the measured frequency of the 
G Raman line of graphite which is 1582 cm -1 (the corre- 
sponding infra red mode is 1586 cm -1 ). This discrepancy 
between Raman and IXS measurements in graphite is so 
far unexplained. 
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FIG. 4: (Color online) Calculated graphene phonon disper- 
sion from DFT (lines) vs. IXS measurements on graphite from 
[U (filled dots), H3l (triangles), and [H (open dots). The 
highest optical branch near K is "corrected" to include GW 
effects following [27l . |28| . and is plotted with a thicker gray 
(red) line. The dashed line is the same branch calculated from 
standard DFT, without GW correction. The cross at T is the 
measured Raman G line frequency in graphite (1582 cm -1 ) 



C. Electron-phonon scattering 



computed for second and more distant neighbors. This 
approximation reproduces very well the k and q depen- 
dence of the electron-phonon matrix elements for elec- 
tronic states with k near K and for optical phonons with 
q near T or near K. This was already verified in [32| by 
direct comparison with DFT calculations. 

We define the average square 
between it and it* at K as (.Dp) f for the E2 9 phonon at 
r. (D^)f is the analogous quantity for the phonon 
at K. From Eqs. IB41 IB5I from App. IB 11 after some al- 
gebra, (D r ) F = 9/4( m ) 2 and (D& F = 9/2( m ) 2 ( m 
is defined in the previous paragraph and the notation 
is consistent with 27]). It follows that, within TB, 
{Dj S }f/{Dj,)f = 2 (that is, this ratio does not depend 
on the actual value of the TB parameter 771). This last 
relation is well reproduced by DFT calculations, within 
LDA or GGA, but not by GW ones (see Table I of 
27]). As a consequence, a single value for 7^ could 
be used to describe reasonably well the DFT electron- 
phonon interaction for phonons in all the Brillouin zone. 
On the contrary, we need two distinct values for r)i, 
rj\ = 5.25 eV/A and 77^ = 6.55 eVA, to reproduce the 
GW value of (Djj-)f and (Dp)jf, respectively, from Table 
I of [27| . Here we will use 771 = 77J 7 for phonons near T 
(those associated to the D' and 2D' lines), and 771 = 7;^ 
for phonons near K (D, 2D, and D + D"). A change of 
r)f and 77^ values will affect the present calculations as 
an uniform intensity scaling of some peaks with respect 
to others. 



D. Electron-light scattering 

Explicit expressions for the D in and D out matrix el- 
ements are given in App. IB 31 We assume that the po- 
larization of the incoming and scattered light are on the 
graphene (x,y) plane. The computed Raman intensity 
Ii.o depends on two indexes determined by the polar- 
ization of the incident (i = x,y) and of the scattered 
light (o = x, y). The polarizations are chosen so as to 
reproduce different kind of Raman experiments. In the 
parallel polarization case, the incident and scattered light 
are parallel polarized and I\\ = I xx + I yy . In the trans- 
verse polarization case, the incident and scattered light 
are perpendicularly polarized and ij_ = I xy + I yx . If 
the light is not polarized I unp oi = Ixx + I yy + Ixy + Iyx- 
Unless specified differently calculations are done in the 
non-polarized case. In Sec. IIII B 31 the effects of parallel 
and transverse light polarizations are discussed. 



The electron-phonon scattering matrix elements 
AiJ q _„ are obtained from TB (explicit expressions are 
given in App. IB 2p and depend on the parameter 771, de- 
fined as the derivative of the nearest-neighbors hopping 
parameter with respect to the bond length. The present 
approach neglects the derivative of the hopping param- 
eters (with respect to the atomic positions) for hopping 



E. Electron-defect scattering 

Defect scattering is treated within the Born approx- 
imation. Namely, the defect scattering operator Hp is 
the difference between the TB Hamiltonian in presence 
of the defect and that of the defect free system. Ho is 
determined by considering three distinct kind of defects. 
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i) The on-site defects: defects that change the value of 
the on-site TB parameter by SV . 

ii) The hopping defects: change the value of one of the 
first-neighbor hopping TB parameters by St\. 

hi) The Coulomb defects: charged impurities adsorbed 
at a distance h from the graphene sheet that interact 
with graphene with a Coulomb potential. Following [33j , 
we consider an environment dielectric constant k = 2.5. 

We remark that these are very simplified prototypical 
models and that a realistic description of a given type 
of impurity, which is beyond the present scope, will re- 
sult in a combination of these three kind of perturba- 
tions. However, it is reasonable to expect that the present 
three models describe the most important characteristics 
of certain kind of defects. For instance, the on-site de- 
fect is the most simple description of an hydrogen atom 
bound to a carbon atom in the graphene sheet. Hopping 
defects are any defects that lead to deformations of the 
carbon-carbon bonds in graphene. A Coulomb defect 
describes any charged atom or molecule adsorbed over 
the graphene sheet. Explicit expressions of the three de- 
fect scattering operators Hd are given in App. IB 41 The 
three models are characterized by the parameters SVq, 
Sti, and h, whose values will be specified in the discus- 
sion. The results will be expressed as a function of the 
defect concentration rid — N^/Aq, where A Q = \f2>/2af ) 
is the graphene unit-cell area, being ao — 2.46 A the 
graphene lattice spacing. 

Note that the Raman intensity of the defect-induced 
lines (e.g. D, D', and D") is proportional to the aver- 
age number of defects in the unit cell, Nd (Eq. [3]) . This 
is because the scattering from defects on different sites 
is considered as incoherent, which is reasonable for low 
defect-concentrations. In particular, for on-site and hop- 
ping defects, the defect-induced intensities are propor- 
tional to a on = nd(SVo) and to Uhopp — nd(dti) , being 
nd the defect concentration. Through the text, we will 
specify the value of these parameters, in order to make 
meaningful the comparison of the defect-induced line in- 
tensities with those of the phonon-phonon lines (e.g. 2D, 
2D', and D' + D"). 



We consider y as the sum of two contributions 

a a(ep) a(D) 

7k =7 k +7 k • 



(7) 



The first is due to electron-phonon scattering. It is an 
intrinsic broadening (present in perfectly crystalline sam- 
ples) and, according to the Golden rule, is 



oc(ep) 

7 k 



2ir 



(8) 



where a refers to tt or tt* bands, the sum is performed on 
a uniform grid of N q q points in the Brillouin zone and 
on all the phonon branches v. A good approximation 
of 7 a ( ep ) is obtained by considering conic bands (|e| = 
hvpk, being Vp the Fermi velocity) and only the two 
phonons E2 9 at T and A^ at K, with energies fiwr and 
fi qjK- By defin ing (g*) = y/h/(2Mwr)(Dl) F and (g& = 
^h/(2MLu K )(D^) F (see Sec. EE]), Eq. M becomes: 



7, 



a(ep) 



[2<<?r)^a(|e| - far) + (gl)N a (\e\ - fa K ) 



«■<«>- £(i£:)W 



(9) 



where N a is the electronic density of states of the a = 
tt or tt* band, being ao the lattice spacing and 8{x) the 
Heaviside step function. Using the parameters of the 
present work, N a (e) = O.O79O8eV~ 2 |e|0(|e|) and for |e| > 
0.196 eV 



7±2 = 41 -89( 



e\ - 0.1645) meV, 



(10) 



where e is expressed in eV. 

The second contribution in Eq. [7] is due to electron- 
defect elastic scattering. It is extrinsic (it is induced by 
the presence of impurities and depends on the sample 
quality) and is 



a(D) 

7 k 



N *W £ l< k '< a ^ k ' a ^ 5 ^ - £ £'), (11) 



F. Electronic linewidth 

An electronic state |ka) (a — tt* or 71") has a finite life- 
time (which is associated to a line broadening energy 
7^ = H/t£) because the electronic states interact, e.g., 
with phonons and with defects. The broadening energies 
7k in the denominators of the DR amplitudes K ( e.g. in 
Eqs. |31 E]) are the sum of the broadenings of the corre- 
sponding electronic states. As examples, in both Eqs. |H 

El 7k = 7if + 7k: 7k 3 . = 7k+q + 7£, and 7 £ = 7l f + 7 £. 
For a — 7r* or tt, 7^ is the full-width at half maximum 
of the electron/hole spectral function as measured, e.g., 
by ARPES. 



where the sum is performed on a uniform grid of N' k k' 
points in the Brillouin zone. The electron-defect scatter- 
ing operator Hd is defined as in Sec. Ill El and App. IB 41 
and depends on the considered kind of defect. Nd is the 
average number of defects in the unit cell. 

Fig. [5] shows ^ epS > and 7^' for the three kind of de- 
fects we considered (7^ = ■y (on \ = y( ho PP\ or 
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(Coul) 



). The 7 in Fig. [5] are calculated with 
Eqs. [51 [TT] and are plotted as a function of the energy 
of the corresponding electronic state (e£ or e£). y( e P> 
is compared with the conic-band results of Eq. [TO] As 
expected, the two results are similar for energies smaller 
than 1 eV. 

j(° n ) anc i ry(Coui) are uruV ocally determined by the en- 
ergy and, in Fig. are represented by lines. y( on \ in 
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particular, is proportional to the density of states. On 
the contrary, j^ ep ' and ^ Ho pp"> display a dispersion asso- 
ciated to the fact that different k electronic states with 
the same energy can have a different life-time. However, 
the dispersion is relatively small, and for the present pur- 
pose they will also be considered a function of the energy. 
All the contributions (7 (ep) , 7 (on) , ^ }io pp"> and 7 ( cw )) 
increase with energy and display a noticeable asymme- 
try between positive and negative energies due to the 
graphene electron/hole asymmetry. 




Energy (eV) 

FIG. 5: (Color online) Electronic linewidth as a function of 
energy, (a) Contribution of electron-phonon scattering to 
the electronic linewidth, ^ ev \ compared to conical bands 
results (Eq. imp , (b) Contribution of on-site and hopping 
impurity scattering to the electronic linewidth. 7'°™' is pro- 
portional to a on — nd(SVo) 2 and ^ hovv ^> is proportional to 
cthopp = nd(Sti) 2 (Sec. Ill E|) . We, thus, plot 7^' /at, where the 
label "i" refers to "on" (on-site defect) or to "hopp" (hopping 
defect), (c) Contribution of Coulomb impurity scattering to 
the electronic linewidth. ricoui is the Coulomb impurity con- 
centration. The distance between graphene and the charged 
impurity h = 0.27 nm (see the discussion in Sec. IIII D"| . 

In actual calculations (e.g. in Eqs. [H [5]) we neglect 
the dependence on k and we use 

i£ = iE = i£ = i tot , (12) 

where 7*°* depends only on the excitation energy e^, on 
the kind of defect D and on its concentration no, through 

f = ^)( e[ )+f)( £i ,„ D ). (13) 

7 are the sum of the two contributions for n an ir* bands 
in a small energy range close to half the excitation energy 
e L . As an example, 7^ = 7 (ep) (ez/2) + ^ ep \-e L /2), 
where ^ ep \e) is the average of ^ ep ^ from Fig. [5] at that 
energy, in particular, for cl > 1-0 eV, 

7 (ep) (e L ) = ( 18.88 e L + 6.802 e| ) meV, (14) 

where £l is expressed in eV. While comparing these val- 
ues with literature, notice that j( tot ) and the 7's corre- 
spond to the sum of the width of electrons and holes and 



are, thus, roughly two times bigger that the width of elec- 
tronic states. To give some examples, for = 2.4 eV, 
and for the typical defect concentrations of the present 
work, a on = a hopp = 6.4 x 10 13 eV 2 cm~ 2 , 7 (on) = 5 meV 
and j( h °PP) = 12 meV, and for n Co ui = 10 12 cm" 2 , 
^(Coui) = 01 meV. On the other hand, for e L = 2.4 eV, 
^(ep) _ me y j g tfi e dominant contribution and, in 
several cases, we will just consider j tot ~ j( e P>. Similar- 
values of 7* ot ~ ^( e P) have been extracted from measure- 
ments in [34[ (note that je-ph of [34| corresponds to 
7^ ep ^/4 in the present notation). 

Finally, in charged graphene a further contribution to 
the broadening due to electron-electron interaction (34J 
can be relevant when 0.06 1 | > 7(^/4 where 6f is the 
Fermi energy (see e.g. Eq.8 of [3J]). For electron/hole 
concentrations of the order of 10 12 cm -2 this contribution 
is negligible and, here, it is not considered. 

III. RESULTS AND DISCUSSION 

This section presents the calculation of the double res- 
onant (DR) Raman spectra of graphene and discuss the 
results. Sec. IIII Al describes the overall agreement with 
measurements. Sec. IIII Bl describes the dependence of 
the spectra on excitation energy and light polarization. 
Sec. IIII CI describes the dependence of the Raman in- 
tensities on various parameters such as the electronic 
linewidth, the excitation energy, and the defect concen- 
tration. Sec. IIII Dl describes the dependence of the spec- 
tra on the type of defect. Sec. IIII El is dedicated to the 
interpretation of the results. It is focused on some spe- 
cific issues such as the determination of the most relevant 
processes and phonons, the role of quantum interference, 
and on the interpretation of the small width of the main 
DR Raman lines. 



A. Overall agreement with measurements 

Figs. [5] and [7] compare the present calculations with 
Raman spectra of Refs. [H, HH, for an excitation energy 
€l = 2.4 eV. In Fig. [SJ below 2000 cm -1 the processes 
are due to phonon-defect scattering and calculations are 
done considering only the hopping defects (this choice 
is justified in Sec. IIII Dp . using the parameter a^opp — 
6.4 x 10 13 eV 2 cm~ 2 (see Sec. Ill E[ ), which reproduces the 
measured ratio of the integrated areas between D and 
2D lines of [11] . Above 2000 cm -1 , all the processes are 
due to two-phonon scattering. We remark that the G line 
is a single-resonant process which is not included in the 
present calculations. 

The agreement between calculations and measure- 
ments is extremely good. In particular, all the lines 
observed experimentally, even the small intensity ones, 
are present in the calculated spectra and the relative in- 
tensities among phonon-defect lines (such as the D and 
the D') or among two-phonon lines (such as 2D, 2D' , or 
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FIG. 6: (Color online) Intensity vs. Raman shift for el 
— 2.4 eV. Comparison of the present calculations with the 
measurements from Notice that our model includes only 
double- resonant processes and, thus, the G line is not present. 
Measurements correspond to a defect concentration — 10 12 
cm' 2 . Calculations are done using j tot — 96 meV, and hop- 
ping defects with cthopp = 6.4 x 10 13 eV 2 cm -2 . All the inten- 
sities are normalized to the maximum value of the 2D line. 




FIG. 7: (Color online) Intensity vs. Raman shift for el — 
2.4 eV. Comparison of the present calculations with the mea- 
surements from Ref. 4!]. The figure reports only two-phonon 
processes. Calculations are done using r ) tot = 84 meV. All the 
intensities are normalized to the maximum value of the 27? 
line. The inset shows the D+D" band in a different scale. 



As far as the line frequencies are concerned, calcula- 
tions and measurements display some small deviations of 
the order of a few meV. We remark that the line frequen- 
cies are determined by a subtle interplay between the 
phononic and electronic energy dispersions, and that the 
present dispersions are obtained from state of the art ab- 
initio computational methods which correctly reproduce 
ARPES and IXS measurements (Sec. Ill 51) . A correction 
of the electronic or of the phononic dispersions, to repro- 
duce with more precision the Raman frequencies, would 
be done at the expense of introducing fitting parameters 
to the model, which is beyond the present scope. 

B. Dependence of the spectra on the laser 

This section describes the dependence of the spectra on 
excitation energy and light polarization. Excitation en- 
ergies vary from 1.2 to 4.0 eV, which are energies mainly 
used in actual experiments. 

1. Dependence of the main lines on the excitation energy 



4 



3 

>> 

a; 2 
3 









• 


i 1 - 
1 




D 


D' 


e L = 1.2eV 


2D' 


.D 


2D 




Jl 




e L =2.4eV D+E) .. J 


l 2D' 
^ 




D 






'I 


l. 




2D 








E L =3.8eV D± p" | J 


L, . 2 P' 



1500 2000 2500 3000 

Raman Shift (cm ) 



FIG. 8: Calculated Raman spectra for el = 1.2 eV and r y tot 
= 32 meV, e L = 2.4 eV and j tot = 84 meV, e l =3.8 eV 
and 7 <ot = 170 meV. Calculations are done using hopping 
defects with othopp = 6.4 x 10 13 eV 2 cm -2 . All the intensities 
are normalized to the corresponding 2D line maxima. 



D + D") are correctly reproduced. The most remarkable 
agreement relates to the line widths. Indeed, the present 
model reproduces very well the measured small widths 
of the D, D', 2D and 2D' lines. Moreover, the model 
reproduces quite well the symmetric Lorentzian shapes 
of the 2D and 2D' lines and the asymmetric shape of 
D + D" band. We remark that, in the present model, the 
only parameter used to fit the Raman data is tXhopp- This 
parameter determines the ratio of the D vs. 2D inten- 
sities but does not affect the relative intensities among 
phonon-defect or among two-phonon lines, the width of 
the lines, and their shape. 



Fig. |8] displays the calculated spectra of the main dou- 
ble resonant Raman lines for three different excitation 
energies. In all cases, we use the electronic broadening 
rytot _ j(ep) ^ ca i cu i a ted at the corresponding excitation 
energy (Sec. Ill F|) . In general, by increasing the excita- 
tion energy, the bands become broader and the relative 
intensities change. The behavior of the 2D line is partic- 
ularly interesting. At er, = 2.4 eV, the 2D line presents 
a Lorentzian lineshape with a relatively small linewidth, 
while at er, — 3.8 eV, it is much broader showing two 
components with smaller, 2D~ , and higher, 2D + Raman 
shifts, as discussed in detail in Sec. IIII E 41 Here, we just 
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FIG. 9: (Color online) Phonon dispersion of graphene along 
high symmetry lines. Bold crosses indicate the phonons that 
mostly contribute to the D, D' , D" , D 3 , D 4 and D 5 Raman 
bands, for e_L = 2.4 eV. Dotted crosses indicate phonons that 
also contribute to the D, D' , D 3 , and D 4 bands but with 
smaller intensity. The crosses are determined from the maxi- 
mum of Iq as defined in Sec. IIII E 21 
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FIG. 10: (Color online) Raman shift as a function of excita- 
tion energy. Upper panel: two-phonon bands. Lower panel: 
disorder induced bands. Our results compared to experimen- 
tal data from Ref. |9| (circles) and Ref. [35j] (triangles). 



remark that the presence of a small width 2D line with 
Lorentzian shape is commonly used to detect a graphene 
monolayer in samples containing flakes with a different 
number of graphene layers [l[ . According to Fig. [H this 
kind of experiment makes sense only when it is done at 
(-L 2.4 eV, but not at higher excitation energies. 

Fig. [UJ shows the wavevector and the branch of the high 
symmetry phonons which mostly contribute to the DR 
graphene lines, for el = 2.4 eV. The figure display the 
phonons associated with the single-phonon Raman lines 
D, D', D", D 3 , D 4 and D 5 , where D 3 , D 4 and D 5 re- 
fer to the small intensity lines of Fig. [TT] The D line is 
associated to the phonon branch affected by the Kohn 
anomaly (thick grey line in Fig. [4]). This branch, near T, 
becomes almost transverse (TO). The D' line is associ- 
ated to the branch which, near T, is almost longitudinal 
(LO). The two-phonon bands, such as the 2D, 2D' and 



D + D" are associated with the emission of two phonons 
which, in the scale of Fig. [3J are almost indistinguishable 
from those of the D, D', and D" lines. 

Fig. [TU] shows the calculated shift of the main Raman 
lines as a function of the excitation energy, ej,. The Ra- 
man shift of the D and 2D lines increases with increasing 
laser energy. The D' Raman shift does not show a mono- 
tonic behavior but it does not change significantly. The 
D + D" Raman shift is almost constant for el between 
1.2 and 1.8 eV, and decreases for e L >1.8 eV. Fig. [TU] 
also shows the experimental data from Ref.Q for the 
2D and D + D" lines and from Kef. for the 2D line. 
The good agreement with measurements is not surprising 
since the dispersion of a DR line as a function of cl is 
determined by the phonon dispersion and in Ref. [271 it 
was already shown that the present phonon dispersions 
(obtained from DFT plus GW corrections) reproduce the 
measured D line shift as a function of ez,. The behavior 
of the shift as a function of ej, is easily understood by 
comparing with the phonon dispersions in Fig. |UJ For 
instance, for the D line, when the excitation energy in- 
creases, the phonons mostly involved in the DR process 
move away from K, and their frequencies are higher. The 
same reasoning explains the behavior of the D' frequency. 
For the two-phonon lines, one has to consider the fre- 
quencies of the two phonon involved. For instance, the 
2D line Raman shifts are twice as large as the D ones. 
For the D + D" line, the energy of one phonon branch 
increases, while the other decreases while moving away 
from K. 



2. Small intensity bands 
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FIG. 11: Calculated Raman spectra for small intensity bands. 
Calculations are done using el = 2.0 eV and 7 °* = 65 meV 
(upper), e L = 2.4 eV and 7*°' = 84 meV (middle), e L = 2.8 
eV and 7 = 106 meV (lower). We consider hopping defects 
with a^opp — 6.4 x 10 13 eV 2 cm -2 . All the intensities are 
normalized to the corresponding 2D line maxima. 

The calculated spectra display some small intensity 
bands which are shown in Fig. [TTJ Some of these bands 
are extremely weak and it is not clear whether they could 



11 



2400 



2200 



« 2000 
c 

ts 

| 1800 

"g 1000 

* 800 

f 600 
53 

I 400 
200 



1 

2D" 




— - 




i 


D+D5 






— * 






ir+p4_ 


D'+D3 , 


* r- 


• 

• — 


I 








i 


1 1 1 ' 1 ' 1 




U 

D5 












D4 












m 




i.i.i.i 



2.0 2.2 2.4 2.6 

Excitation energy (eV) 



FIG. 12: (Color online) Raman shift vs. excitation energy for 
the small intensity bands of Fig. 1111 Upper and lower panels 
display results for two-phonon and defect-induced bands, re- 
spectively. Upper panel calculations arc compared with mea- 
surements from [3q ] (dots) and [37| (diamonds). 



be possibly measured, on the other hand the D" is ob- 
served @, [H| and the bands that we label as D' + D 4 and 
D' + D 3 have been measured recently [H,|33|- Fig.[T2"lrc- 
ports the shift of these small intensity bands as a function 
of the excitation energy. The agreement with available 
measurements is good. Fig. [5] reports the high symmetry 
phonons associated with the bands that we label as D 3 , 
D 4 , D 5 , and D" . The D 3 and D 4 bands are associated 
with phonons near T, that have a momentum very simi- 
lar to the momentum of the phonons associated to the D' 
line. The D 5 and D" bands are associated with phonons 
near K, with a momentum very similar to the momen- 
tum of the D phonons. The D 3 , D 4 , D 5 , and D" bands 
are however much weaker than the D and D' ones, be- 
cause the electron-phonon coupling (between ir electronic 
bands) for those branches, is much weaker than the one 
of the D and D' (see pf). 
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FIG. 13: (Color online) Comparison of calculated Raman 
spectra done with different light polarizations. Calculations 
are done using el = 2.4 eV and j tot = 84 meV (upper plot), 
or el = 3.8 eV and -y tot = 170 meV (lower plot). We used 
hopping defects with a^opp = 6.4 x 10 13 eV 2 cm -2 . The in- 
tensities are normalized to the corresponding 2D line maxima 
calculated with unpolarized light. "Parallel" and "transverse" 
refer to In and I± as defined in Sec. Ill Dl 



3. Dependence on the light polarization 

So far, we have shown calculations done with unpolar- 
ized light. We now discuss how the results are affected 
by the use of polarized light. For parallel and trans- 
verse polarizations, we calculated I\\ and I± as defined in 
Sec. Ill Dl Fig. Q2] compares the results obtained for cl 
= 2.4 eV and £l = 3.8 eV. The intensity in the parallel 
polarization case is considerably larger than in the trans- 
verse one, as expected. For el = 2.4 eV, the spectrum 
shape almost does not depend on the polarization and the 
ratio I\\/I±_ is about 2.7, in reasonable agreement with 
measurements in graphite [HI, graphene [3^| and earlier 
theoretical predictions For ez — 3.8 eV, the D and 
2D bands split into two components (see Sec. IIIIE"4l for 
a detailed discussion) and the intensity ratio between the 
two components depends on the polarization. For exam- 



ple, the intensities of the two components of the 2D band, 
2D + and 2D~ , are very similar within transverse polar- 
ization, while the 2D + intensity is slightly higher than 
the 2D~ one, within parallel polarization. This finding 
is very remarkable since it could lead to measurable ef- 
fects. 



C. Dependence of the Raman intensities on the 
various parameters 

In this section we discuss how the intensity of the main 
DR Raman lines is affected by the various parameters 
such as the electronic linewidth (Sec. IIII C ip . the exci- 
tation energy (Sec. IIII C 2[) . and the defect concentration 
(Sec. IIII C 3[) . In general, the absolute value of the inten- 
sities is affected by these parameters, however, we will 
mainly focus on how the ratio of the intensities of dif- 
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ferent lines is affected, since this last quantity can be 
measured more easily. 

1. Dependence on the electronic broadening 
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FIG. 14: (Color online) Integrated areas under the 2D, 2D' 
and D + D" lines [A{2D), A{2D'), and A(D + £>")] as a 
function of the electron + hole linewidth (7*°*), for el = 
2.4 eV. The areas are normalized to A(2D) calculated with 
j tot = 7 (ep) = 84 meV. For clarity, A(2D') and A(D + D") 
are multiplied by 20. Symbols are calculations, lines are the 
fit from Eq.Hni Inset: A(2D)/A(2D') ratio. 

As already discussed in Sec. Ill Fl the broadening pa- 
rameter 7 tot (the sum of the electron and hole linewidths, 
see Eq. IT2|) results from an intrinsic component (due to 
electron-phonon scattering), which depends on the laser 
energy, and from an extrinsic component which increases 
by increasing the defect concentration. Eventually, in 
charged (doped) graphene, a further contribution due to 
electron-electron scattering can be relevant. The actual 
value of 7 tot , which depends on the defect concentration, 
determines in a measurable way also the intensities of the 
two-phonon lines (which are not defect induced). Indeed, 
Fig. Q3] reports the integrated areas under the 2D, 2D' 
and D + D" lines [A(2D), A(2D'), and A(D + D'% as a 
function of 7*°*. The areas of these lines decrease by in- 
creasing 7*°*. In general, for all Raman lines studied here, 
the intensity decreases when the electronic linewidth in- 
creases, at fixed defect concentration. This is because, in 
Eq. [TJ an increase of the imaginary values ij tends to kill 
the double resonance condition. 

It is interesting to notice that also the ratio of the two 
areas, A(2D)/A(2D'), depends on 7*°* (inset of Fig. HI). 
This result is particularly appealing since the ratio of the 
two areas can be measured in a relatively easy way. The 
measured value of A(2D)/A(2D') compared to the inset 
of Fig. UM (which is obtained for ej, = 2.4 eV) could, thus, 
be used to determine experimentally the electron+hole 
linewidth 7*°* and, in particular, its components due to 
defects and/or to electron-electron scattering in doped 
samples (keeping in mind that for large doping the value 



of the electron-phonon interaction itself is expected to 
change [4(| and, thus, the inset of Fig. UM cannot be used 
as it is). For j tot — j( ep > = 84 meV, which is suitable 
for comparison with pristine graphene, A(2D)/A(2D') 
= 21.5, in agreement with experimental works which re- 
ported A(2D)/A(2D') as being 27 and 26 ± 3 HJ. 

In [23J it has been shown that, if the electronic bands 
can be considered conic, the dependence of A(2D) and 
A(2£>') on 7*°* should be A= A /(7*°*) 2 , where A is 
a constant. This functional form, however, cannot be 
used for a quantitative description of the present results. 
Indeed, the integrated areas as a function of 7*°* reported 
in Fig. [14] can be fitted by a similar, but different, law: 

A(2D) = 9374/(( 7 tot ) 2 + 48.5 2 ) 
A{2D') = 629/(( 7 tot ) 2 + 80.0 2 ) 
A(D + D") = 438/(( 7 tot ) 2 + 59.6 2 ), (15) 

where 7*°* is expressed in meV. An explanation of the dis- 
crepancy between Eqs.[T5land the model of [23[ (which is 
based on a simplified description of the electronic bands) 
is probably associated to the importance of a proper in- 
clusion of the trigonal warping and of the electron/hole 
asymmetry in the description of the electronic bands 
(Sec. HH}. Another result of [H| is that 

A(2£>)/A(2£>') = 2(^/^)4 x {u2D ,/ U2D f_ ( i 6) 

Eq. [16] is obtained by rewriting the equation in the last 
paragraph of [23| using the notation of Sec. Ill CI and 
considering lu2d and u>2D' are the frequencies associated 
with the two Raman lines. Indeed, for large 7 tot , the 
ratio A(2D)/A(2D') from Eqs. [TBI does not depend on 
7* ot . However, using the parameters of the present work, 
Eq. EH gives A(2D)/A(2D')=6.8 which is almost two 
times smaller than A(2D)/A(2D / )=14.7 obtained from 
the limit 7'°* — !• 00 of Eqs. [15] This second discrepancy 
with the model of [23] is so far unexplained, since in 
this limit the effect of electron-hole asymmetry should 
become negligible. We also remark that the model of 
[H predicts that the ratio A(2D)/A(2D') does not de- 
pend on the excitation energy e^. In the following we 
will show that, on the contrary, A(2D) / A(2D') strongly 
depends on ej,. 

2. Dependence on the excitation energy 

The intensity of the 2D line decreases by increasing 
the excitation energy (Fig. [T5|). The most important 
contribution to the decrease comes from the fact that the 
electron/hole broadening 7'°' increases by increasing ex,. 
This can be deduced from Fig. [T5] which also shows the 
results for a fictitious system in which 7*°* is kept to a 
fixed value independent from ex- Indeed, in this second 
case, the dependence of A(2D) on ex is much less marked 
than in the full calculation. 

Fig. [ToT a) reports the calculated ratio of the inte- 
grated areas under the bands, A(2D')/A(2D) and A(D + 
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FIG. 15: (Color online) Integrated area under the 2D line as a 
function of the excitation energy el- The defect concentration 
is zero. The full line is obtained by including the dependence 
of the broadening on e L , 7*°* = 7 {ep] '(ex,) (see Sec-HTF). The 
dashed line is from an unrealistic simulation in which j has 
been kept fixed to a constant value 7 iot — ^( ep ) 
84 meV, independent from el- 
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FIG. 16: (Color online) Ratio of the integrated areas under 
Raman bands as a function of excitation energy, (a) Two- 
phonon bands: our results compared to experimental data 
from Ref. [l). (b) Disorder induced bands from hopping im- 
purities, with cthopp = 6.4 x 10 13 eV 2 cm~ 2 . 



D")/A(2D), as a function of the excitation energy ex. 
These ratios considerably change in the range of exci- 
tation energies of the figure. A(2D')/A(2D) decreases 
and A(D + D")/A(2D) increases rapidly. The values 
calculated for e L = 2.4 eV compare reasonably well with 
those obtained from the measurements of [lj. In the 
last paragraph of Sec. IIIIC II we discussed the model of 
|23j , which was used to theoretically determine the ratio 
A(2D')/A(2D). The simplified model of [H predicts 
that the ratio A(2D')/A(2D) does not depend on el- 
On the contrary, from Fig. I16f a) , this dependence is very 
important. Using Eq. \Tfi\ (which is adapted from [23j]) 
and using, for consistency, the parameters of the present 
work, one obtains A(2D')/A(2D) = 0.15. This value is 



significantly higher than 0.09, which we obtain for the 
smallest e L of Fig. ITBTa). 

Fig. Iil)r b) reports the ratio of the integrated ar- 
eas under the defect-induced bands, A(D') / A(D) and 
A(D")/A(D). Here, we consider again only hopping im- 
purities. We also remark that the present approach is 
expected to be valid in the limit of small defect concen- 
tration. For small excitation energies the D" band inten- 
sity is very small in comparison to the D one. For larger 
excitation energies the D" relative intensity increases, 
reaching A(D")/A(D) = 0.09 when e L = 4.0 eV. On the 
other hand, the intensity of the D 1 band compared to 
the D band decreases by increasing the excitation en- 
ergy. For el up to about 3.0 eV the D 1 band is more 
intense than the D" band, while for el > 3.2 eV, the D" 
is slightly more intense than the D', 
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FIG. 17: (Color online) Intensity of the D and 2D Raman 
lines as a function of the defect concentration for el = 2.4 eV. 
Calculations are done using hopping defects and are reported 
as a function of the parameter cthopp = nd(Sti) 2 {rid is the 
defect concentration and 5t\ the hopping parameter), in the 
upper horizontal scale. The lower horizontal scale is obtained 
by considering Hi = 8.0 eV. (a) Id is the maximum of the 
intensity of the D line; symbols are experimental data from 
0| • The dashed line is a linear fit of the Id calculated values 
for Ud < 5 x 10 11 cm -2 . Theoretical and experimental in- 
tensities have been normalized by their maximum values, (b) 
Integrated areas under D and 2D bands, A(D) and A(2D). 
Experimental data are from [Tl|. Theoretical and experi- 
mental areas are normalized by A(2D) at minimum defect 
concentration. The vertical line indicates the defect concen- 
tration of 7 x 10 12 cm" 2 (a hopp = 4.5 x 10 14 cm" 2 eV 2 ) for 



which the two contributions to the electronic broadening are 
equal: j^ ' 



We now discuss how the intensities of the Raman bands 
are affected by defect concentration rid- We recall that 
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two-phonons Raman lines (such as the 2D) depend on 
rid only through the electronic broadening parameter 
jtot (jZq. [12]) . -f tot is given by the sum of an intrin- 
sic component y( ep > (due to the electron-phonon inter- 
action) and an extrinsic defect-induced component y( D ' 
which increases linearly by increasing rid (Eq. 1 13[) - On 
the other hand, the defect-induced Raman lines (such as 
the D line) depend on rid through two distinct mecha- 
nisms. First, it depends on rid through 7 * as for the 
two-phonon lines. Second, there is a proportionality fac- 
tor between the Raman intensity and the the number 
of defects in the sample (/ oc Nd in Eq. [3]). Basically, 
for a higher number of defects there are more scattering 
events that can activate the defect-induce lines, which, 
in crystalline samples, are not Raman active. In the fol- 
lowing discussion, we will consider only hopping defects. 
As already shown in Sects. Ill El and III Fl the calculated 
Raman spectra depend on the defect concentration, rid, 
only through the parameter ahopp — n-d(8ti) 2 , being 5t\ 
the hopping parameter. 

Fig. Unreports the D line peak maximum (Id) and the 
integrated areas under the calculated 2D and D lines, 
A(D) and A(2D), as a function of the parameter ahopp, 
for e L = 2.4 eV. For a hopp = 4.5 x 10 14 cm- 2 cV 2 , the two 
contributions to the broadening are equal, 7^ = ^i e p) _ 
The corresponding ahopp is indicated in Fig. [T7] with a 
vertical line. The intensity of the 2D line (which corre- 
sponds to a two-phonon process) monotonously decreases 
by increasing the defect concentration. For small defect 
concentrations (ahopp < 10 14 cm~ 2 eV 2 ) 7^ <C 7^), 
ry tot j( e P) slightly depends on the defect concentra- 
tion, and A(2D) is almost constant. For higher defect 
concentrations, 7'-°) becomes the dominant contribution 
to 7*°*, which, as a consequence, becomes more sensi- 
tive to the defect concentration. The increase of 7*°* by 
increasing the defect concentration is associated to a de- 
crease of A(2_D), because of the mechanism discussed in 
Sec. IIIIC II 

The intensity of the D line (which is a defect in- 
duced process) has a different behavior. For low de- 
fect concentrations, it increases almost linearly, then it 
reaches a maximum, and finally decreases. This behav- 
ior results from the interplay of two competing mech- 
anisms. For small defect concentration 7^' <C ^ ep ^ 
and 7* ot ~ j( e P) . I n this region, the intensity is ex- 
pected to increase linearly (/ cx Nd in Eq. Indeed, 
the calculated intensity is well reproduced by a linear 
fit up to ahopp < 10 14 cm _2 eV 2 (compare the continu- 
ous line with the dashed one in Fig. [TT1 upper panel). 
For ahopp > 4.5 x 10 14 cm~ 2 eV 2 , the dependence of the 
broadening 7*°* on the defect concentration becomes the 
dominant mechanism, leading to a decrease of the inten- 
sity as for the 2D line. It is remarkable that the defect 
concentration for which ahopp = 4.5 x 10 14 cm _2 eV 2 (ver- 
tical line in Fig. I17|) almost coincides with the maximum 
value reached by the D intensity, Id. 

Fig. [T7] compares calculations with the intensities of 
the D and 2D measured in 0, [H| as a function of the 



defect concentration. So far, we have discussed theo- 
retical results as a function of ahopp — nd(Sti) 2 . ahopp 
defines the upper horizontal scale in Fig. [T7J To make 
the comparison with measurements we need to attribute 
a value to the hopping energy 8t\. The best fit to mea- 
surements is obtained for 8t\ = 8.0 eV. This value is used 
only to rescale the horizontal axis of Fig. [T7] and defines 
the defect concentration as reported in the lower horizon- 
tal axis of Fig. [TTJ The measured behavior as a function 
of the defect concentration is well reproduce by calcula- 
tions. It is remarkable that the same value Sti = 8.0 eV 
can be used to fit equally well the D and the 2D line 
data. The value Sti = 8.0 eV is very high. However, one 
should notice that in Ref. @, E2 defects were induced in 
graphene by means of Ar + ion bombardment. This tech- 
nique leads to the formation of Carbon multi-vacancies 
in the sample. In Ref. I?}, the defect average size is es- 
timated, by means of scanning tunnel microscopy, to be 
1.85 nm. On the contrary, the present model considers 
only point defects (the hopping parameters is changed by 
Sti for a single isolated carbon-carbon bond) . The large 
value 5ti = 8.0 eV is, thus, to be considered as an effec- 
tive variation of the hopping parameter that mimics the 
existence of an extended defect (a realistic description of 
the defect should be done by considering the variation 
of the hopping parameters associated to many different 
neighboring sites). For less damaging defects, Sti will be 
smaller and the critical defect concentration, above which 
the D line intensity begins to decrease, will be larger than 
that of Fig. [T7] 

Finally, the behavior of the D line intensity as a func- 
tion of the defect concentration has been discussed in 
literature using different models @, EH (see also [42l|). 
To make a comparison, it can be useful to restate the 
present finding as follows. According to the DR pertur- 
bative model, the intensity of the defect-induced lines 
decreases by increasing the defect concentration when 
7^ becomes higher than y( e P\ that is when the average 
length an electron/hole travels in between two scatter- 
ings events with a defect becomes smaller than the aver- 
age length an electron/hole travels before scattering with 
an optical phonon. 



D. Dependence of the spectra on the type of defect 

Here, we discuss how the results depend on the type 
of defect. Calculations were done using three different 
model defects namely, hopping defects, on-site defects, 
and Coulomb ones (see Sec. Ill El for a description of the 
relevant parameters). Fig. [18] compares calculations with 
the measurements from [111 ], which correspond to a de- 
fect concentration rid — 10 12 cm~ 2 and 6l = 2.4 eV. For 
the hopping and on-site defects, the calculations are done 
using ahopp — a on — 6.4 x 10 13 eV 2 cm -2 , which, for the 
hopping defect, reproduces the ratio between the inte- 
grated areas of the measured D and 2D lines of [ll[ . By 
choosing Sti = SVq — 8.0 eV (see also the discussion in 
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FIG. 18: (Color online) Calculated Raman spectra obtained 
for three different kind of defects (hopping, on-site, and 
Coulomb), compared with the measurements of [ll[ done 
at e_L = 2.4 eV. The Raman G line is not described by the 
present model. Calculations are done using j tot — 96 meV. 
Other relevant parameters are given in the text. All intensi- 
ties are normalized by the corresponding 2D maximum. The 
intensity of the Coulomb impurity spectrum is enhanced by 
10 2 for clarity. 



Scc. lIII C~3j) , the above values of a correspond to a defect 
concentration rid — 10 12 cm~ 2 . For Coulomb impurities, 
the distance between the impurity and graphene is h = 
0.27 nm and rid = 10 12 cm -2 . 

From Fig. [TH the hopping defect is the best model 
to study defect-induced Raman processes. Indeed, con- 
trary to the other models, the hopping defect provides 
a ratio of the intensities of the D and D' lines which 
is in good agreement with measurements. The intensity 
ratio between D and D' strongly depends on the kind 
of model defect, suggesting that this ratio could pos- 
sibly be used to experimentally determine the kind of 
defects present in a graphene sample. From Fig. 1181 we 
also notice that Coulomb defects (charged impurities out- 
side the graphene plane) provide an almost undetectable 
contribution to the Raman signal. Indeed, for a defect 
concentration of rid = 10 12 cm~ 2 , the D line is absent 
and the D' intensity is almost three orders of magni- 
tude smaller than the experimental one. We recall that 
Coulomb defects could be an important source of scat- 
tering during electronic transport in graphene (see (43j 
and refs. therein). The fact that they are not detectable 
by Raman spectroscopy (which is routinely used to char- 
acterize experimentally the quality of graphene samples) 
is, thus, a relevant issue which deserves some more com- 
ments. 

The present simulations consider a very short 
graphene/impurity distance h, in order to enhance the 
Raman signal of the Coulomb impurities. Indeed, h = 
0.27 nm is the distance between K atoms and graphene 
planes in the KCs intercalated graphite. This distance 
corresponds to the experimental conditions of [43} , where 



K + ions are deposited on graphene. In the case, where 
the impurities are charges trapped in the substrate (e.g. 
SiC-2) a longer distance (e.g. 1 nm) is more appropriate. 
It is not surprising that the contribution of Coulomb im- 
purities to the D line is completely negligible. Indeed, 
the Fourier transform of the Coulomb potential is maxi- 
mum close to r and decays as 1/q far from it, Eq. IB101 
and the D line is due to phonons near to the K point 
and far from T. This argument, also, suggests that the 
D' band, which is due to phonons near T, should be more 
sensitive to the presence of Coulomb impurities. Accord- 
ing to calculations, this is actually the case, however for 
cl = 2.4 eV and rid = 10 12 cm~ 2 the ratio of the inte- 
grated area A(D')/A(2D) = 1.5 x 10~ 4 , meaning that 
the presence of a D' band due to Coulomb impurities 
should not be detectable. The use of smaller energy laser 
increases the intensity of the D' signal since the excited 
phonons are nearer to T. However, for cl = 1.2 eV and 
n d = 10 12 cm" 2 , A(D')/A(2D) = 8.0 x 10~ 4 , which is 
still very small. Within the present model, A(D') / A(2D) 
increases linearly by increasing the impurity concentra- 
tion, rid- rid, however, cannot be higher than 10 14 cm~ 2 , 
which corresponds the density of K atoms in KCs- On the 
other hand, for Coulomb impurity concentrations higher 
than 10 12 cm" 2 doping effects should become important. 
These should be associated to an increase of the electron- 
electron scattering contribution to the electronic broad- 
ening 34], which, in turn, will prevent the D' intensity to 
become detectable. Concluding, the presence of charged 
impurities is not associated to a Raman D band. A D' 
band is present, but should not be easily detectable. 



E. Interpretation of the results 

This section is dedicated to the interpretation of the 
results. Sec. IIIIE II describes which are the most impor- 
tant processes associated to the DR. Sec. MI E"2l describes 
which are the phonon wavevectors contributing to each 
Raman band. Sec. 1111 E 31 analyzes the dominant direc- 
tions of the phonon wavevectors and Sec. IIIIE 41 is dedi- 
cated to the interpretation of the small width of the main 
DR Raman lines. 



1. Dominant Processes and Interference Effects 

In this section we analyze which are the dominant pro- 
cesses among those described in Fig. [TJ We distinguish 
between two classes of processes: processes aa are those 
in which the two intermediate scattering processes are 
associated to both electron states or to both hole states 
(namely the processes eel, ee2, hhl, and hh2, using the 
notation of Fig. [T]); processes ab are those in which the 
two scattering processes are associated one to an elec- 
tron state and the other to a hole state (ehl, eh2, hel, 
and he2 in Fig. [I}. The distinction between aa and ab 
processes holds for both phonon-defect and two-phonon 
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lines. 

In general, for all the simulations performed here, the 
ab processes are, by far, dominant over the aa ones, that 
is, the largest part of the Raman intensities are due to ab 
processes. This is true for both phonon-defect and two- 
phonon lines. In general, among the ab processes, all the 
four processes ehl, eh2, hel, and eh2 are associated to in- 
tensities of the same order of magnitude. Indeed, Fig. [Hfl 
shows a typical Raman spectrum, in which we compare 
the actual spectrum Itot with two spectra obtained by 
including only aa processes, I aa , or ab ones, I O0 . More 
precisely, Itot is the Raman intensity computed includ- 
ing all the processes; I aa is computed by restricting the 
sums in a and /3 in Eqs.[3]only to the eel, ee2, hhl, and 
hh2 processes; l a b is computed by restricting the sums 
in a and /3 in Eqs. [3]only to the ehl, eh2, hel, and he2 
processes. In general, I tot ^ l aa + l a b- From Fig. [191 
lab ~> laa for both the D and the 2D lines. 

The dominance of the ab processes is due to quantum 
interference effects. In particular, from Eq.[31 the Raman 
intensity for a given q results from a sum over k of K (k) 
scattering amplitudes, which are complex numbers. The 
sum of these complex numbers can interfere in a con- 
structive way, as for the ab processes, or in a destruc- 
tive way, as for the aa processes. In particular, the DR 
condition determines that for some resonant electronic 
wavevectors k r , |if(k r )| should have a maximum. This 
maximum can be enhanced or suppressed by the inter- 
ference of K(k r ) with the K (k) at wavevectors k which 
are not exactly at the resonance (this point is further 
discussed in App. [D]). It is important to remark that, 
according to the present calculations, the DR scattering 
amplitudes K are complex numbers in which the real and 
imaginary parts are of the same order of magnitude even 
for the k = k r wavevectors that satisfy the DR condition. 

To quantify the importance of quantum interference, 
we consider a fictitious Raman intensity I, which is ob- 
tained by substituting their modulus \K\ to the scat- 
tering amplitudes K in Eqs. [3] As example, in Eqs. [3] 



we substitute JSP 



E M ^(k.i) /N k , with \pp 



Ek,^l^(k,q)| /N k . 

Thus, within the intensities I, the presence of possible 
destructive interference effect is cancelled. Fig. Q1J] shows 
a typical I spectrum, in which we compare l aa and I ao 
obtained by solely including aa or ab processes. The ra- 
tio l a b /laa is very different from I ao /I a a f° r both the D 
and the 2D lines. In particular, I ao is no more dominant 
and it is always comparable in intensity to l aa . Thus, 
the fact that l a b l aa is indeed due to destructive inter- 
ference effects. Moreover, certain lines of the fictitious I 
spectrum, such as the D 1 or the 2D', do not appear as 
narrow and well defined lines as they are in the actual 
Raman spectrum, I. Thus, interference effects also play 
a role in determining the shape of certain lines. 

Notice that, often, when discussing the DR processes, 
it is used a simplified argument which consists in finding 
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FIG. 19: (Color online) The upper panels compare the calcu- 
lated Raman spectrum hot with spectra determined consid- 
ering only aa processes, l aa , or ab processes, l ao . More pre- 
cisely, Itot is determined considering all the processes shown 
in Fig. [1] laa is computed by considering only eel, ee2, hhl, 
and hh2 processes; l a b is computed by considering only ehl, 
ehl, hel, and he2 processes (see the text). The lower panels 
display fictitious Raman intensities 1 obtained by substituting 
to the DR scattering amplitudes K in Eqs. |3j their modulus 
\K\ (see the text). The two lines l aa and l a t are obtained by 
considering only aa and ab processes, as before. Calculations 
are done using el = 2.4 eV, j = 84 meV, and hopping 
defects with a^opp = 6.4 x 10 13 eV 2 cm -2 . All the intensities 
are normalized to the 2D line maximum of hot- 



the electronic and phonon states which let two (or more) 
of the denominators in Eq. [T]go to zero. The assumption 
is that the physics is lead only by those scattering ampli- 
tudes K which satisfy the DR condition. This simplified 
approach, which we call the "resonance argument" , has 
been extensively used in literature with success (e.g. to 
determine the momenta of the phonons associated to cer- 
tain lines), despite the fact that, within this approach, 
the possible role of quantum interference is completely 
neglected. The results of the previous paragraph show 
that in certain specific situation the "resonance argu- 
ment" can be very misleading. For example, on the basis 
of a "resonance argument" one would deduce that the 
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intensity associated aa processes are of the same order of 
magnitude than that associated to the ab ones (indeed, 
iaa ~ lab in Fig. [19]), which is not the case. 



We remark that several authors describe the DR by 
simply consider the aa processes (usually the ee pro- 
cesses in Fig. [U 0]) , as it is done in the seminal work by 
Thomsen and Reich [loj . However, following the present 
conclusions, these processes cannot be used alone to de- 
scribe quantitatively the intensities of the D and 2D 
lines. The importance of interference effects in deter- 
mining the shape of the DR Raman lines has been al- 
ready outlined by Maultzsh et al. in [29j. However, 
Ref. [29[ just consider ee processes and completely ne- 
glects the ab ones, which are the most important. The 
fact the ab processes should be dominant for the 2D line 
has been argued by Basko in Ref. (23|. But, this conclu- 
sion is reached on the basis of a "resonance argument". 
Indeed, according to Ref. [23[ , the ab processes should be 
dominant because within these process one can reach a 
condition in which all the transitions are real (non vir- 
tual) and the three denominators of Eq.Q]can be nullified 
simultaneously (triple resonance). As already said, this 
kind of arguments cannot be applied to describe the in- 
tensity of the 2D line (basically, the conclusion is good 
but the argument is wrong) . The best way to understand 
this point is to put to zero the phonon energies hu) p h in 
all the denominators of the Raman scattering amplitudes 
K (e.g. in Eqs.[HE]). By doing this, the triple resonance 
condition of Basko applies also to the aa processes (not 
only to the ab). However, actual calculations show that 
l a b remains much larger than I oa even when fvjj p h = 0. 
Actually, the intensity and the shape of the 2D line are 
marginally affected by including or not haj p h in the de- 
nominators of the Ks (see Fig. [27] in App. [Cj. We also 
remark that the triple resonance argument does not ex- 



2. Phonons wavevectors associated to the Raman lines 



plain why l ao ^> l aa also for the D line. Finally, Ref. [44 1 
argues that quantum interference in real space plays a 
crucial role in enhancing the role of the ab processes ver- 
sus the aa ones, for the D line. However, the model of 
Ref. [13], predicts a behavior which is in contrast with 
the present calculations [45j . Notice that the model of 
44J was developed to describe extended defects such as 
edges, while here we are considering point defects. 



The main conclusion of this section is that the ab pro- 
cesses (ehl, eh2, hel and eh2 processes of Figs. [TJ SJ) are 
responsible for most of the Raman intensity because of 
quantum interference. We remark that this conclusion 
is not due to the complex details of the present calcula- 
tions but can be deduced with a very simplified model in 
which the scattering matrix elements in the numerator of 
Eq[T]are constant, the phonon energies in the denomina- 
tors (e.g. hujq in Eqs. H] [5]) are neglected, and in which 
the electronic bands are conic. This simple model can 
also be used to shed light on the role played by quantum 
interference, see App. [DJ 
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FIG. 20: (Color online) Decomposition the intensity of the 
most important Raman bands into its components associated 
to phonons with a given wavevector q, X q . The rhombi are 
the graphene first Brillouin zone. For each band, we con- 
sider the contribution to the Raman intensity in a window of 
frequencies corresponding to that particular band 46]. The 
intensities are normalized to the maximum of each band. Cal- 
culations are done using el = 2.4 eV, 7*°' = 84 meV, and 
hopping defects with cthopp = 6.4 x 10 13 eV 2 cm -2 . 



We now discuss which phonons are responsible for the 
lines presented in Figs. |6j and [7] In Fig. [20J we con- 
sider the most important Raman lines and we decom- 
pose the Raman intensity of a given band into its com- 
ponents associated to phonons with a given wavevector 
q. For the defect-induced bands, D, D' and D" , we plot 
Z q = Yltlqi> an d f° r the two-phonon bands, 2D, 2D' 
and D + D", we plot X q = £* If v ^ with I** and J™, 
defined in Eq. 3 and the symbol * indicates that the sum- 
mation is restricted to a frequency window corresponding 
to a given Raman band (see [4q). The q-dependent in- 
tensity Z q discloses which are the phonon wavevectors q 
that mostly contribute to a given Raman line. The most 
remarkable result from Fig. [501 is that these phonons be- 
long to limited regions of the BZ consisting in very narrow 
(almost one-dimensional) lines. As expected, the D, D", 
2D and D + D" Raman bands originate from phonon q 
wavevectors belonging to a closed line around the K and 
K' high symmetry points. 

In literature, the DR condition on the virtual tran- 
sitions is often used to determine the Raman- dominant 
phonon- wavevectors (see, e.g., [3, H, E3, Hi, S3, El ) ■ To 
verify the validity of such a procedure, we focus on the 
2D line, which is mostly due to eh processes fSec. HIIET]) 
and consider an excitation energy = 2.4 eV. The DR 
consists in three processes of excitation, phonon scatter- 
ing, and recombination. The k vectors of the electronic 
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FIG. 21: (Color online) Electron and phonon states relevant 
for the 2D line. The rhombi are the graphene Brillouin zone, 
a) The triangularly distorted contour around K is obtained 
from e£ — e£ = 2.4 eV and represents the electronic states 
near K that are excited by a laser with energy 6l = 2.4 eV. 
The contour around K' = 2K is obtained from e£ — e£ = 
2.06 eV and represents the electronic states near K' that are 
deexcited by the emission of a quantum of light with energy 
el — 2uj p h eV, with uj p h = 1354 cm -1 (half the energy of the 
2D line for et = 2.4 eV). b) q n is one of the vectors such that 
the contour near K translated by q n is tangent to the contour 
near K'. c) I q decomposition of the 2D intensity (same figure 
as the 2D panel in Fig. I20|) . The dashed closed line is defined 
by the ensemble of the q„ vectors, d) The dashed line is the 
same as in c). The thick grey (red) line is the phonon iso- 
energy contour obtained from lo^ = 1354 cm -1 . The relevant 
phonon branch, thick grey line in Fig. [4] is disentangled form 
the other branches as in Fig. 2 of Ref. [23]. Notice that the iso- 
energy contours of electron states (panels a, b) and phonons 
(panel c, d) have opposite trigonal warpings. Notice also that 
phonon iso-energy contours in Fig. 2 of Ref. [28l ] are plotted 
with respect to the K' of the present notation. 




^ Sj — 2hd) ph 



FIG. 22: (Color online) Scheme of the double resonant pro- 
cess associated to the 2D line. The momenta of the phonons 
mostly involved are indicated as "inner" and "outer". 
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FIG. 23: (Color online) Lower panel: momenta of the inner 
and outer high symmetry phonons which mostly contribute to 
the 2D band. The lines are obtained from the vectors connect- 
ing the isoenergy electronic contours corresponding to that 
excitation energy. For example, the values for cl = 2.4 eV 
are the moduli of the "inner" and "outer" vectors reported in 
Fig [21] The symbols are obtained from the maximum inten- 
sity in the X q plots (as those in Fig. [20] or in the left panels of 
Fig.l26[l corresponding to that excitation energy. Upper panel: 
frequency of the "inner" and "outer" phonons reported in the 
lower panel. 



states which are excited by a laser with energy 6l form a 
triangularly-distorted closed line, as the iso-energy con- 
tour surrounding the K point in Fig. 121a . The states in- 
volved in the emission of a quantum of light with energy 
£l — 2friOq (recombination) form a second triangularly- 
distorted closed line, as the iso-energy contour surround- 
ing the K' point in Fig. [2Tb . These iso-energy contours 
are expected to give the important contribution to the 
DR, although the energy is not conserved in the inter- 
mediate virtual transitions. The intermediate DR pro- 
cesses are associated to a phonon q and the important 



processes are expected to be those associated to q vec- 
tors that connect the two triangles of Fig. [2Tb . In par- 
ticular, let us translate the K triangle by q and let us 
consider the nesting vectors (q n ) for which the K tri- 
angle becomes tangent to the K' one, as in Fig. [2Tb . 
These phonon- wavevectors are expected to dominate the 
Raman spectra, since for such nesting vectors there is 
a high density of electronic transitions satisfying the DR 
mechanism [l4],|48j]. The q n vectors are shown in Fig. [2Tb 
as a dashed white line which is compared with the Ra- 
man intensity X q from our most precise calculation (as 
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in Fig. |20|) . Within the scale of the figure, the nest- 
ing vectors reproduce very well the maximum of the Z q , 
meaning that the simple picture of Fig. [2Tb provides a 
quantitative prediction of the relevant phonon momenta. 

To generalize the analysis to an arbitrary laser 
excitation-energy, we now consider, the isoenergy elec- 
tronic contours as those of Fig. [2Th for different values 
of €l- For each cl, we determined the phonon q„ vec- 
tors that are nesting the corresponding contours. Among 
these points, we consider only the vectors along high sym- 
metry lines. In this case the nesting vectors, qi nnC r and 
qoutcr, can be easily extracted from the one-dimensional 
electronic-band dispersion along the high symmetry line, 
as schematically shown in Fig. [22] In the lower panel of 
Fig. [23] we report qinnor and q ou tcr obtained by the DR 
condition of Fig. [22] as a function of e^. In Fig. [231 wc 
also report the corresponding vectors obtained by find- 
ing the maximum intensity in the Z q plots (as those in 
Fig. [2TJI) corresponding to that excitation energy. The 
sets of q vectors obtained with these two different proce- 
dures nicely coincide. 

We remark that the simplified scheme of Figs. [2Tb 
and [22] is used for the 2D line, and that its validity 
comes "a posteriori" after the comparison with our most 
precise calculations. The analogous construction for the 
2D' line works equally well, as can be seen in Fig. [2Tb. by 
comparing the nesting vector profile (dashed line) with 
the I q decomposition of the 2D' intensity. 

3. Dominant directions of the Raman phonon-wavevectors 

A close look at Fig. [20] reveals that the most intense 
contributions of D, D", 2D and D + D" are due to q 
points along the high symmetry directions K— >T and 
K'— >r. The D' and 2D' bands originate from a closed 
line around T and the most intense contributions are due 
to q points along the high symmetry r~>M direction. 

To analyze the results we consider the following def- 
initions. The intensities I q (Fig. [2D]) form, basically, 
a closed profile surrounding one high symmetry point 
(K for the D and 2D lines, and T for the 2D'). Tak- 
ing the high symmetry point as the reference, we con- 
sider how the intensity of a given Raman band varies 
as a function of the direction q of the vector q. Thus, 
in the lower panel of Fig. [25] we plot I q = J Q 9 qdqlq, 
where the integral is done in a region containing the 
most intense contribution. It is also interesting to con- 
sider the intensity weighted average phonon frequency 
associated to a given Raman band and to a given q 
point, (w q ). As example, for the two-phonon lines 

K) = ^«+<)]/E^ where the sum - 

mation is restricted to the corresponding frequency win- 
dow [IfjJ]. This quantity, basically, gives the frequency of 
the phonons associated to that Raman band. In analogy 
to I q , we define (uiq) as the average of (u; q ) along a di- 
rection q of the vector q. Here, also, the origin of q is 
K for the D and 2D lines, and T for the 2D'. Fig. [2S] 




FIG. 24: (Color online) Electron and phonon states rele- 
vant for the 2D' line, a) The triangularly distorted con- 
tours around K are obtained from e£ — e£ = 2.4 eV and 
e k — e k — 2.0 eV. They represents the electronic states that 
are excited by a laser with energy el = 2.4 eV and those that 
are deexcited by the emission of a quantum of light with en- 
ergy e_L — 2oj p h eV, with io v h = 1602 cm -1 (half the energy of 
the 2D' line for = 2.4 eV). q n is one of the vectors such 
that the excited-states contour translated by q n is tangent to 
the second contour. The analogous construction arond K' is 
also shown, b) I q decomposition of the 2D' intensity (same 
figure as the 2D' panel in Fig. I20[) . The two dashed (green) 
closed lines (almost indistinguishable in the scale of the fig- 
ure) are defined by the ensemble of the nesting q n vectors 
among K states or among K' states. 



shows the angular dependence of the averaged phonon 
frequency (ojq) for the D, 2D, and 2D' lines (actually, 
the shifts in the upper panel of Fig. [25] are obtained 
after an average on a small angle interval from 8 — A8 to 
6 + A6). 

Let us consider the D and 2D bands. From Fig. [251 
the phonons along the K— >T directions (in literature 
these are usually called "inner" phonons, Fig. I22j) pro- 
vide a contribution which is almost four times higher 
than the one from the K— s-M ones ("outer" phonons). 
Contrary to the present findings, in literature it is usu- 
ally assumed [2, [l(J Ell that the phonons which mostly 
contribute to the D and 2D lines are outer phonons 
(along K— >M). Only very recently some authors have 
outlined the possible importance of the inner phonons 
(K— yT) \UM22\ . The present finding is counter-intuitive 
and stems from the complex behavior of the scatter- 
ing matrix elements in the numerators of Eq. [TJ To 
understand this point, in Fig. [25] we show the results 
of calculations in which the numerators in Eq. [TJ arc 
taken as a constant (that is, independent form k and 
q as, e.g., in Eqs.[H[S]). Within this simplified approach 
(which completely neglects, for example, the dependence 
electron-phonon scattering matrix elements on q) the 
outer phonons become dominant (in Fig. 1251 lower panel, 
the intensity has the maximum along the K— >M direc- 
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FIG. 25: (Color online) Angular dependence of the intensity 
(lower panels) and of the weighted average Raman shift (up- 
per panels) for the D, 2D and 2D' bands. The angles are 
measured taking the horizontal direction in Fig. [20] as ref- 
erence. Thus, for the D and 2D bands, zero degrees is the 
K— !>r direction in the BZ, while ±60 degrees are the K— s-M 
one. For the 2D' band, zero degrees is the T— >K direction, 
while, ±30 degrees are the T— >M direction. In the lower pan- 
els, the solid lines correspond to our most precise calculation. 
Dashed lines correspond to an approximated simulations in 
which the electron-light, electron-phonon, and electron-defect 
scattering matrix elements are kept constant (see the text). 
Calculations are done using e_L = 2.4 eV, y tot = 84 meV, and 
hopping defects with aho PP = 6.4 x 10 13 eV 2 cm~ 2 . 



tion for both D and 2D), in agreement with the simpli- 
fied models previously used in literature, but in disagree- 
ment with our most precise calculations. Concluding, 
inner processes are dominant for both D and 2D lines. 
A proper description of the electronic scattering matrix 
elements (in particular of the electron-phonon coupling) 
is crucial to obtain this result. 



4- The width of the Raman bands 

One of the most interesting feature of the simulated 
Raman spectra of Figs. [6] and [7] is the narrow width 
of the bands, which reproduces the measured spectra. 
The narrow width of the D and 2D lines is indeed sur- 
prising since already at ej, = 2.4 eV the electronic states 
involved in the Raman process display an important trig- 
onal warping (i.e. the electron isoenergy contour are tri- 
angularly distorted as in Fig. [2Tk). In the presence of 
trigonal warping one should expect the excited phonons 
to have energies distributed in a broad range. Indeed, 
previous calculations [3, [TtJ did not reproduce narrow 
lineshape of the DR lines. The present improved descrip- 
tion of the electronic scattering matrix elements partially 
explains such narrow lines. The most important role is 
played by the phonon energy dispersions. The upper pan- 
els of Fig. [55] show that, for the D, 2D and 2D' lines 
at €l = 2.4 eV, the excited phonons have almost the 
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FIG. 26: (Color online) Calculated 2D line for the excitation 
e L = 3.8 eV and y tot = 170 meV. Top right panel: intensity 
vs. Raman shift. The line appears as a broad band with two 
maxima near 2790 cm" 1 (2D") and 2840 cm" 1 (2D+). Left 
panels: mapping of the Raman intensity in the Brillouin zone 
(as in Fig. I20[) of the two components 2D~ 2D + obtained 
by integrating in the corresponding frequency windows [46l |. 
Central panels: angular dependence of the weighted average 
Raman shift and of the intensity, as in Fig. 1251 



same energy (within ^5 cm" 1 ), despite the strong elec- 
tron trigonal warping. This fact explains the small width 
of the DR Raman lines and it is due to the details of the 
phonon dispersion we used. Indeed, with a reasonable 
description of the electronic trigonal warping and using 
a rough description of the phonon energies, larger dis- 
persions in frequencies (and broader Raman lines) are 
found [l7|. Ref. [HI has clearly demonstrated that the 
phonon trigonal warping is important and that it is op- 
posite to the electronic one. The present results show 
that, as already argued in Ref. (28[, the interplay between 
the electronic and phononic trigonal warping provides a 
sort of cancellation. This results in the small dispersion 
of the phonon frequencies of the upper panel of Fig. [231 
and, consequently, in the small width of the associated 
Raman lines. 

To illustrate the concept of trigonal warpings cancella- 
tion, Fig.[2TH compares the line of the nesting vectors q n 
(white dashed line, see Fig. [2Tb and Scc. lIIIE"2|) with the 
iso-energy contour of the phonons having half the energy 
of the 2D (thick red lines). The two lines nicely resemble 
each other, meaning that all the nesting phonons have 
nearly the same energy and, as a consequence, the 2D 
line width is small. If the phonon isoenergy contour was 
different, the two lines would not superimpose and the 
2D line would have a broader shape. The perfect cancel- 
lation of electronic and phononic trigonal warping breaks 
down for laser energy in the UV range. Indeed in the up- 
per panel of Fig. [UJ we report, as a function of cl, the 
frequency associated with the inner and outer phonons. 
At el — 2.4 eV, the frequencies associated to inner and 
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outer phonons almost coincide. However, for higher el, 
the two frequencies become different, meaning that for 
a sufficiently high the 2D line is expected to become 
broader. 

Indeed, according to our most precise calculations, at 
el = 3.8 eV the Raman 2D band appears much broader 
than the one at cl — 2.4 eV and displays two maxima at 
2790 cm" 1 and 2840 cm" 1 (Fig. At e L = 3.8 eV 
(Fig. the angular dependence of the average fre- 
quency shift is more dispersive than in the cl = 2.4 eV 
case. The inner phonons correspond to the highest fre- 
quency components, 2D + at ^2840 cm -1 , and the outer 
phonons to the lowest one, 2D~ at ^2790 cm^ 1 . In 
Fig. [21] we also show the q vectors decomposition of the 
intensities of the 2D + and 2D" components. For the 
2D + , the shape is triangularly distorted and the max- 
imum corresponds to the inner phonons, while for the 
2D~ the maximum corresponds to the outer phonons. 



IV. CONCLUSIONS 

We calculated the double resonant Raman spectrum 
of graphene with a computational method which tries to 
overcome the most common approximations used in liter- 
ature. Calculations are done using the standard approach 
based on the golden rule generalized to the fourth per- 
turbative order [lCf (Eq. [1]). We determined the Raman 
lines associated to both phonon-defect processes (defect- 
induced excitations of q^O phonons, such as in the D, 
D' , and D" Raman lines) and two-phonons processes (ex- 
citations in a defect-free sample of a -q and a q phonons, 
such as in the 2D, 2D', or D + D" lines). The lowest- 
order processes (excitation of a q=0 phonon, such in 
the G line) and higher-order processes (such as in the 
D + D' line at ~2900 cm -1 , which is usually attributed 
to a defect-induced excitation of two phonons q and q' 
with q+q'7^0) are not described by the present approach. 

The electronic summation is performed all over the two 
dimensional Brillouin zone and all the possible phonons 
(with any wavevector) are considered. Electronic bands 
are obtained from a 5-neighbors tight binding (TB) ap- 
proach in which the parameters are fitted to reproduce 
ab-initio calculations based on density functional theory 
(DFT) corrected with GW. This procedure provides a 
Fermi velocity (the slope of the Dirac cone) in good agree- 
ment with measurements and a good description of the 
trigonal warping. The resulting electron/hole asymmetry 
is not negligible. The phonon dispersion is obtained from 
fully ab-initio DFT calculations corrected with GW. This 
procedure is necessary to obtain a good description of the 
slope of the phonon branch associated with the D and 
2D lines, near K. The electron-phonon, electron-light, 
and electron-defect scattering matrix elements are ob- 
tained within the TB approach. The defect-induced Ra- 
man processes are simulated by considering three differ- 
ent kinds of model defects: i) on-site defects, obtained by 
changing the on-site TB parameter; ii) hopping defects, 



obtained by changing one of the first-neighbors hopping 
TB parameters; iii) Coulomb defects, corresponding to 
charged impurities adsorbed at a given distance from the 
graphene sheet, which interact with graphene through a 
Coulomb potential. 

The electronic linewidth (the inverse of the electronic 
lifetime), which turns out to be a very relevant parame- 
ter, is calculated explicitly considering the contributions 
from electron-phonon and electron-impurity scattering. 
To give an idea, for = 2.4 eV, in the absence of de- 
fects and for zero doping, the sum of the electron and 
hole linewidths is 7*°* = 84 meV (which is roughly two 
times the FWHM of the electron spectral function). 

By looking at the overall shape of the typical Ra- 
man spectra, for an excitation energy of cl = 2.4 eV, 
the agreement between calculations and measurements is 
very good. In particular, all the Raman lines observed ex- 
perimentally, even the small intensity ones, are present in 
the calculated spectra and the relative intensities among 
two-phonon lines (such as 2D, 2D', or D + D" lines) or 
among phonon-defect lines (such as the D and the D' 
lines) are correctly reproduced (being the hopping de- 
fect the best model to study defect-induced Raman pro- 
cesses). The most remarkable agreement between the- 
ory and measurements relates to the line widths. In- 
deed, the present calculations reproduce very well the 
measured small widths of the D, D' , 2D and 2D' lines. 
Moreover, calculations reproduce quite well the symmet- 
ric Lorentzian shapes of the 2D and 2D' lines and the 
asymmetric shape of D + D" band. We remark that, in 
the present model, the only parameter used to fit Ra- 
man measurements, cthopp, determines the ratio of the D 
vs. 2D intensities but does not affect the relative inten- 
sities among phonon-defect or among two-phonon lines, 
the width of the lines, and their shape. 

We determined how the Raman spectra change by 
changing the laser excitation energy from 1.2 to 4.0 
eV, which are the energies mainly used experimentally. 
All the visible lines change in position, intensity and 
shape. In particular, the 2D line has a small-width 
Lorentzian shape for cl < 2.4 eV and it is asymmetric 
and broader at = 3.8 eV. The measured shift of the 
Raman line position as a function of cl is well reproduced 
for all the available measurements. The calculated spec- 
tra also display some small intensity bands associated to 
acoustic phonons. Some of them, such as the D' + D 3 and 
the D' + D 4 (in the 1800, 2000 cm" 1 range) are actually 
visible in the measured spectra [36|, [37| • Finally, for high 
energy excitations, e.g. 6l — 3.8 eV, the most intense 
Raman lines (2D and D) change shape and intensity as 
a function of the polarization of the light. This finding 
is remarkable since it could lead to measurable effects. 

We determined how the intensity of the main DR Ra- 
man lines is affected by various parameters such as the 
electronic linewidth, the excitation energy, and the de- 
fect concentration. The absolute intensity of the double 
resonant Raman lines is strongly affected by the actual 
value of the electronic linewidth, 7'°'. In general, the 
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intensity of a DR Raman line decreases when the elec- 
tronic linewidth increases (at fixed defect concentration) 
because the electronic broadening tends to kill the dou- 
ble resonance condition. According to the present find- 
ings, also the ratio of the intensities of the 2D and 2D' 
lines depends on 7*°*. This result is particularly appeal- 
ing since the measurement of this ratio (followed by the 
comparison with the present calculations) could be used 
to determine experimentally the electron/hole linewidth 
7* ot and, in particular, its components due to defects 
and/or to electron-electron scattering in doped samples. 
We determined how the intensity ratio among various Ra- 
man lines change as a function of the excitation energy of 
the laser. In particular, we determined the evolution of 
A(2£>')/A(2£>), A(D + D")/A(2D), A(D")/A(D), and 
A(D')/ A(D) [where A(X) is the integrated area under 
the X Raman line] as a function of the excitation en- 
ergy. All these ratios considerably change in the range 
of excitation energies available experimentally, however 
measurements to compare with are not presently avail- 
able. 

We studied the dependence of the D and 2D lines in- 
tensity on the defect concentration, comparing to recent 
measurements @, HH • We first remind that the electronic 
linewidth 7 ' is given by the sum of an intrinsic compo- 
nent 7"( ep ) (due to the electron-phonon interaction) and 
an extrinsic defect-induced component j( D *> which in- 
creases linearly by increasing the defect concentration. 
The intensity of the 2D line monotonously decreases by 
increasing the defect concentration n^. Indeed, the 2D 
line (which is a two-phonon process) depends on rid only 
through the electronic linewidth j tot i which, in turn, in- 
creases by increasing n ( i- The intensity of the D line has 
a non-monotonic behavior. The D line (which is a de- 
fect induced process) depends on rid through two distinct 
mechanisms: first there is a proportionality factor be- 
tween the Raman intensity and rid, second, the linewidth 
7 tot depends on rid as for the 2D line. For small rid, 
ry tot j( e P) and the D intensity increases linearly with 
rid- For high rid, the dependence of 7*°* on rid becomes 
the dominant mechanism, leading to a decrease of the 
intensity, as for the 2D line. The maximum of the D 
intensity is reached for the defect concentration corre- 
sponding to the condition y( D ) ~ ^( e P). 

We have compared Raman spectra calculated with the 
three different model defects. The intensity ratio between 
the defect-induced D and D' lines strongly depends on 
the kind of model defect, suggesting that this ratio could 
possibly be tuned in actual experiments by selecting spe- 
cial kind of impurities on the sample. Charged impurities 
outside the graphene plane (Coulomb defects) could be 
an important source of scattering during electronic trans- 
port. However, according to the present calculations, 
they should provide an almost undetectable contribution 
to the Raman signal, the D line being completely ab- 
sent and the D' having an intensity orders of magnitude 
smaller than the 2D line. 

Finally, the analysis of the results has focused on cer- 



tain specific issues currently debated. 

Among the different possible DR processes, the 
electron-hole ones (processes in which both electronic and 
hole states are involved in the scattering, ab in the text) 
are responsible for most of the Raman intensity of both 
the D and the 2D lines. Several authors (e.g. [l(|) de- 
scribe the DR by simply considering electron-electron or 
hole-hole processes (processes in which only electrons or 
only holes are involved in the scattering, aa in the text) 
which, according to the present findings, give a negligi- 
ble contribution to the Raman intensity. The dominance 
of the electron-hole processes stems from the presence of 
a destructive quantum interference that kills the contri- 
bution of the electron-electron and hole-hole ones. This 
conclusion is not due to the complex details of the present 
calculations but can be deduced with a very simplified 
model, easy to implement. 

The most intense contribution to both the D and 2D 
lines is due to phonons along the high symmetry direc- 
tions K— > r (inner phonons) . This is contrary to the 
common assumption [l|, 03, [H| that the phonons which 
mostly contribute to the D and 2D lines belong to the 
K— >M direction (outer phonons). The present result 
(the dominance of the inner phonons) is counterintuitive 
and stems from the complex behavior of the electronic 
scattering matrix elements in the numerator of the dou- 
ble resonance scattering amplitude. 

The observed small width of the 2D line at — 2.4 eV 
is explained as a consequence of the interplay between the 
opposite trigonal warpings of the electron and phonon 
dispersions: the excited electronic states form a trian- 
gularly distorted profile having vertex along the K— >M 
direction, while the phonon isoenergy contour is a trian- 
gularly distorted profile having vertex along the K— s- T 
direction. Because of this, the excited phonons (both the 
inner and the outer ones) have almost the same energy 
and, as a consequence, the 2D line-width is small. At 
higher excitation energies this condition is no more ver- 
ified and the 2D line becomes broader and asymmetric. 
For instance at ez, = 3.8 eV the calculated spectrum dis- 
plays two maxima corresponding to a main component at 
^2840 cm" 1 (due to inner phonons) and to a less intense 
one at ^2790 cm -1 (due to outer phonons). 
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Appendix A: Raman double-resonant scattering 
amplitudes 

Explicit expressions are now given for all the dou- 
ble resonant scattering amplitudes K pd (k, q, v) and 
K pp (k, q, is, fj,), which have been included in the sums 
of Eq. [3] The following processes are described diagram- 
matically in Fig. [T] The arguments k, q, v, and [i are 
dropped for simplicity. The sign ± before each K is de- 
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termined by the fermionic statistics of the carriers. The 
broadening energies 7k in the denominators of the DR 
scattering amplitudes K are the sum of the broadenings 
of the corresponding electronic states (see Sec. Ill F|) . As 



examples, in K% el 7^ = 7J 
7k = 7k* +7k- InA7 d ~ A 
7k = 7£+ q + 7£ +q - 



hi 7 k = 7 k 



7£, 7if = 



7k +q H 

= 7Cq" 



There are eight phonon-defect (pd) processes. 
Process eel: the electron is first scattered by a phonon and then by a defect, 
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Process ee2: the electron is first scattered by a defect and then by a phonon, 



K. 



pd 



(k7r|D out |k7r*)(k7r*|Aif q)i/ |k- q) 7r*)(k-q, 7r*|.ff D |k7r*)(k7r*|A„|k^) 



Process hhl: the hole is first scattered by a phonon and then by a defect, 

rpd (k7r|D out |k7r*)(k-q,7r|ff D |k7r)(k7r|Aff q , v |k-q,7r)(k7r*|An|kir) 



K, 



hhl 



{e L - ef +el- hw\ 



)(e L 



;2£ 



){e L 



+ el 



Process hh2: the hole is first scattered by a defect and then by a phonon, 

(kn\D out \kn*) (k + q, tt| AH^kn) (kn\H D \k + q, tt) (k7T*| An|k7r) 



K 



hh2 



(e L 



el - hu> v _ q - i\){e L - + el 



k+ q 



-){e L 



Process e/il: first the electron is scattered by a phonon and then the hole by a defect, 



K 



pd 
eh 1 



(k + q, n\D out \k + q, ir*)(kir\H D \k + q, 7r)(k + q, tt*] Aif q)l/ |k7r*)(k7r*| D in \kir) 



(tL - e£; q + ^ +q - ^- q - i\){zL - e£; q + el - 7^ q - i^)( £i - ef + e£ - i^-) 
Process eh2: first the electron is scattered by a defect and then the hole by a phonon, 

rpd (k - q,7r|£> out |k - q,7r*)(k7r|Aff q ,„|k - q, 7r)(k - q, 7T*|ff g |k7r*)(k^*|An|k7r) 



K 



eh2 



(e L - eZ" _ i- f 



k- q T c k- q _ ^-q - ' 



Process hel: first the hole is scattered by a phonon and then the electron by a defect, 

(k - q, n\D out \k - q, tt*) (k - q, 7r*|ff D |k7r*) (k7r| AH^k - q, tt) (k^* | D m \kn) 



K pd 



fa ~ £ k-q + ^k-q - ^-q 



i^)fa-e£* 



C k-q 



Process /ie2: first the hole is scattered by a defect and then the electron by a phonon, 

(k + q,7rp „t|k + q ,7r*)(k + q ,7r*|AH q>l/ |k7r*)(k7r| J ffD|k + q ,7r)(k 7 r*|An|k7r) 



K p - -- 



[e L 



c k+ q c k+ q 



W7 



)fa 



f k+q 



f2t 



)fa 



There are eight two-phonon (pp) processes. 
Process eel: the electron is first scattered by the -qz/ phonon and then by the qp, one, 



Kill 



(eL 



(k7r| D out \kir*) (kit* I Aff_ q , M |k + q, tt*) (k + q, tt* \ AH^kw*) (kir* \ D ln \kir) 

Cq + e k-^-q- j2 f)(^- £ r+ e k-^ 



Ofa 



7k, 
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Process ee2: the electron is first scattered by the qp phonon and then by the -qv one, 

_ (k7r|J> out |k7r*)(k7r*|Ag q , y |k-q,7r*)(k-q,7r*|Ag_ q ^|k7r*)(k7r*|An|k7r) 

ee2 C b A 

(e L - ef + el - fia£ q -fin*- i\){e L - e£l q + ej - fc^ - i\){e L - e£ + el - 
Process ftftl: the hole is first scattered by the -qv phonon and then by the q/x one, 

(k7T|D out |k7T*) (k - q, 7T| Atf _ q , M |k7T) (k7r| Atf q ,„ |k - q, 7T) (kTT* | An|k7r) 



(e L - ef +el- fe^ q - - if )(e L - e£ + e£_ q - fe^, - if )(e x - e£* + cj - if) ' 
Process ft.ft.2: the hole is first scattered by the qn phonon and then by the -qv one, 



(k7T|£> ouf |k7T*) (k + q, tt\ Ag q ,„ |k7T) (k7r| Aff_ q;A1 |k + q, tt) (kTT* \D in \kir) 

(e L - el' +el- fc^ q - fe^ - if )( ei - e£* + e£ +q - Hu% - if )(e L - e£ + e£ - if) ' 

Process e/il: first the electron is scattered by the -qv phonon and then the hole by the q/x one, 
K 



„„ (k + q, 7:\D ou t\k + q, 7r*)(k^| Aff_ q , M |k + q, tt) (k + q, tt* | Aff q ,„|k7r*) (k^* \D rn \kir) 

(e L - e£; q + e£ +q - fe^ q - ^ - if )(e L - e^ q + ej - fe^ q - if ){e L - ef + e£ - if) 

Process eh2: first the electron is scattered by the qp phonon and then the hole by the -qv one, 



(k - q,7r|£> ottt |k- q,7r*)(k7r|Ag q!i/ |k - q, 7r)(k- q, 7r*|Aff- q ^|k7r*)(k7r*|£> m |k7r) 
- ejl q + e£_ q - fi^ q - foj* - if )(e L - e^ + el - tu^ - if )(e L - e£ + el - if) ' 

Process hel: first the hole is scattered by the -qu phonon and then the electron by the q/z one, 

(k - q, 7r|L> out |k - q, 7r*)(k - q, tt* | Air_ q ,„|k7r*><k7r|A£r qi „|k - q, 7r)(k7r*| AnM 



(e L - e-l q + e£_ q - fe^ q - fo^ - i\){e L - ef + e£_ q - fo^ q - i\){e L - e£* + e£ - i^) 
Process /ie2: first the hole is scattered by the qp phonon and then the electron by the -qv one, 

k pp (k + q, ir\D out \k + g, 7r*)(k + q, 7r*|Ag q , y |k7r*)(k7r|Ag_ q , M |k + q, 7r)(k7T*| An|k7r) 

(e L - e£; q + e£ +q - fe^ q - fo^ - if )(e L - ef + e£ +q - ^ - if )(e L - ef + - if) ' 



Appendix B: The Tight-Binding Model 

Here we describe the tight-binding model which is used 
to calculate the electronic structure, the electron-phonon, 
the electron-light and the electron-defect scattering ma- 
trix elements. 



1. Electronic structure 

Let us call \l,s) the orthonormalized p z orbital of the 
s atom (in graphene s = 1,2), in the position r s , in the 
cell identified by the lattice vectors R/ (I = 1, oo). Let us 
consider the wavefunction (normalized in the unit cell) 

|k, s ) = ]T e * k -( R <+^|z, s >. 
i 

Given the tight-binding Hamiltonian H, -ffk.s.s' = 
(k, s|-ff|k, s')/N (N is the number of cells in the crys- 



tal) is the 2x2 matrix: 
where 

/( k ) = - h e e * kc > - 1 8 e eikc? - ** e eikcf 

i=l,3 i=l,3 i=l,6 

5 (k) = -t 2 E e ikC ' -is E eJkCf - f*(k). (B2) 

i=l,6 i=l,6 

Here, U is the i-th neighbor hopping parameter. are 
the three vectors connecting the s = 1 atom with its three 
nearest neighbors (i = 1,3). More in general, C; are the 
vectors connecting the s = 1 atom with the i-th atom in 
the j-th neighborhood. 
By diagonalizing i/k,s,s', 

E -^k : s,s'«kV = e k a ks> 
s'=l,2 
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one obtains the eigenvalues e£ (a — tt,it*) and the eigen 
wavefunctions |k, a) — J2 S a kJk ; s ) : 

e r=. 9 (k)+i/(k)i , «r = ^(^)) 

e ^ = .g(k)-|/(k)| , a^ = -^M (k) ),(B3) 

where 0(k) = /*(k)/|/(k)|. 

Finally, here the overlap matrix is the identity because 
of the use of orthonormal orbitals. In alternative, 
a precise description of the bands can also be obtained 
by using pristine (non-orthonormal) p z orbital with only 
three neighbors interaction parameters at the expense of 
using a non-diagonal overlap matrix (see e.g. (49l. |50|). 



3. Electron-light scattering 

The electron-light interaction is calculated as 

e L 

where P, n and P oll t are the polarizations of the incident 
and scattered radiation, ViT(k) is the gradient of the 
TB Hamiltonian and is a 2x2 matrix, ez, is the incident 
laser energy and e™* is the scattered radiation energy 
(e£ u * = €l — fiU-q for a if pd (q, z/) process and e£ u * = 
£l — fio;^ q — TttJq 1 for a if pp (q, v, tt) process). 

4. Electron-defect scattering 



2. Electron-phonon scattering 



Given a phonon mode qi^, with pulsation w q „ and po- 
larization (s = 1,2 is an atomic index and c = 
1,3 is a Cartersian coordinate index, e^l is normal- 
ized to 1 in the unit cell, corresponding to a displace- 
ment e «.c e iq (R(+-r 3 ) Q f i-^g g a ^ om m ^j^g i unit-cell), the 

electron-phonon scattering matrix element is 



(k + q,a|A/V|k,£) 



2MWq v 



k+q,k a k' 



(B4) 



where M is the carbon mass. All the unit cells give the 
same contribution and the bra-ket integration is done on 
the unit cell (with this choice the numerators of the scat- 
tering amplitudes are independent from the number of 
cells of the crystal) . The 2x2 matrix AH^ k is the 
derivative of the TB Hamiltonian with respect to a pe- 
riodic displacement (with periodicity q) of the atom s 
along the c Cartesian coordinate. By defining rji as the 
derivative of the nearest-neighbor hopping parameter t\ 
with respect to the bond length, 



A ^+ q ,k 



Mk) 



AH, 



2.c 



k+q,k 

Mk) 



K(k + q) 

-V3i 



Mk + q) 
m I h*(k) 



(B5) 



4 = 1,3 



where C\ c is the Cartesian component along the c direc- 
tion of C?", and ao is the graphene lattice spacing. 



We consider three distinct kind of defects. The 
electron-defect scattering operator is defined accordingly. 

i) The on-site deject changes the on-site TB parameter 
of the atom t\ by SVo, in this case we will use the notation 
Hd = V on and 



(ka|y oi 



|ka) = — . 



(B7) 



a = tt or 7r*. Here we have considered t\ in the origin 
and here the bra-ket integration is done all over the space, 
ii) The hopping deject changes the hopping parameter of 
two nearest-neighbor atoms connected by the vector G\ 
by 5t\. H D = Vhopp and 



(ka\Vhopp\b'a) 



fk') 



(B8) 

In the calculations of the 



where <fi is defined as in Eq. 
Raman scattering probability averages among the three 
different C\ vectors are taken. 

iii) The Coulomb defect is a Coulomb impurity with 
charge e, placed at a distance h from the graphene sheet. 
In this case, Ho = Vcoui- The Coulomb potential in the 
position r in the graphene's plane is 



Vcw(r) 



1 



47re K Vr 2 + h? 47re K 



-kh 



(B9) 

where eo the vacuum permittivity, k an environment di- 
electric constant, and the integral is performed on all the 
reciprocal space. By assuming that the p z orbitals are 
localized with respect to ao and h (this is done to avoid 
the introduction of new parameters in the model), 



(ka\Vcoui\k'a) = 



\ e" 



|k-k'+G|ft 



2e Q KA ^ |k - k' 



G| 



1 + e i ( k - k ' +G ) c i0*(k)^(k')l (BIO) 
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where the sum is done on the reciprocal lattice vectors 
G and A is the unit-cell area. 

Note that in the three cases the Raman intensity is cal- 
culated by Eqs. [5] and [3] As a consequence, for the cases 
of on-site and hopping defects the intensity is propor- 
tional to a on = UciSVq and ahopp — ndStf, respectively, 
being n& the impurity concentration. On the other hand, 
for the Coulomb impurities, the intensity is proportional 
to iid, but it also depends on the impurity-graphene dis- 
tance, h, as in Eq.(BlO) above. 



Appendix C: Role of the phonon energies in the DR 




Raman Shift (cm ) 



2600 2800 3000 3200 
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FIG. 27: (Color online) Comparison of a typical Raman spec- 
trum (full calculation) with a test calculation in which the 
phonon energies in the denominators of the DR scattering 
amplitudes K are considered zero. Calculations are done us- 



ing e L = 2.4 eV, Y 



96 meV, and hopping defects with 



Qhopp = 6.4 x 10 eV cm . All the intensities are normal- 
ized to the 2D line maximum value of the full calculation. 



The Raman spectra depend on the phonon frequencies 
Wq through the energy conservation between the initial 
and the final states (expressed in the S functions in Eq. [2]) 
and through the denominators of the DR scattering am- 
plitudes K (e.g. in Eqs. [4j [5j . We performed a serie of 
test calculations in which we consider the phonon ener- 
gies Wq = in all the denominators of the amplitudes 



K (e.g. 



in Eqs. [H [SJ. It turns out that, 



qualitatively, the Raman spectra are not affected. For 
example, the 2D line intensity is basically unchanged, 
while the D one remains of the same orders of magnitude 
(Fig.[27|l. We also checked that the results of Sec. IIIIE1I 
are not affected by the actual value o f in the denom- 
inators. Using the notation of Sec. IIII E 1| by letting 
cjq = in the K denominators, I a f, 3> I aa and I aa ~ lafc 
for both the 2D and the D lines. That is, the ab processes 
are still, by far, the dominant ones. 

Appendix D: A Simple model 



In Sec. IIII E II we have shown that the largest part of 
the DR Raman spectrum is due to the processes involv- 
ing the scattering of both one electron and one hole (ab 
processes). We now show that the same conclusions are 
reached by considering a simple model in which the scat- 
tering matrix elements in the numerator of EqQ]are con- 
stant, the phonon energies in the denominators (e.g. hio^ 
in Eqs.lHE]) are neglected (see discussion in App.JCj, and 
in which the electronic bands are conic: e£ = ±fti>^|k|, 
where vf is the Fermi velocity and k=0 corresponds to 
the high symmetry K point. 

For a given excitation energy cl, the scattering cross 
section associated to a phonon of momentum q are 
I aa (q, €l) and I a &(q, ej,). As usual, aa refers to the eel, 
ee2, hhl, and hh2 processes, and ab to the ehl, e/i2, hel, 
and he2 ones. By using the equations of App. \X\ one 
obtains, 



Wq, e L ) 
K aa (k,q, e L ) 
K ab (k,c[, e L ) 



d 2 k 

WT 2 



K aa (k, q, e L ) 



Iafe(q,£i) 
1 



rf 2 k 



K ab (k, q, e L ) 



(e^ — 2hvpk — i^) (eL — hvp\k + q| — hvpk — i^) (cl — 2hvpk — i^) 

1 

(e L - 2hv F |k + q| — i^) (e L — hv F \k + q| — hvpk — i^) (e L — 2hvpk - 



(Dl) 



In analogy to Sec. IIII E II l aa and l a b are obtained 
by considering only the modulus of the integrand, e.g. 
I QQ = |J d 2 k/(2n) 2 \K aa \\ 2 . Fig. reports the intensi- 
ties / thus obtained for a fixed value of cl, as a function 
of q (the results do not depend on the direction of q). 
As expected from the DR picture, 1(g) has a maximum 
at q = €L/{hvF)- Even with this simplified model, one 



recover the result that ab processes are by far dominant: 
lab 3> Iaa from Fig. 1281 The importance of quantum in- 
terference effects is understood by considering that the 
intensities I a b and I aa (in which quantum interference 
effects arc artificially canceled, Sec. IIII E 1|) are very dif- 
ferent from I O0 and I aa . In particular, I ao and l aa have 
the same order of magnitude. As already noticed in [29| 
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the shapes of 1(g) and 1(g) are very different, thus the 
fact that 1(g) is associated to a well defined narrow line 
is a direct consequence of quantum interferece. Notice 
that, however, the authors of [2i| consider only the aa 
processes. 



To further explain the concept of quantum interference 
we consider that for a fixed value of ez the resonance 
condition g r = cl/(Hvf) (<? = 2 in Fig. |2"5)) , implies that 
the maximum of the intensities are 



FIG. 28: (Color online) Numerical solution of Eqs. ID ll using 
e L = 2.4 eV, 7 = 84 meV, and hv F = 6.49 eVA. q = 2qhv F /e L 
is an adimensional momentum and q = 2 corresponds to the 
double resonance condition. I aa is magnified by 10 2 for clarity. 
I aa and lab are intensities in which quantum interference has 
been artificially suppressed (see the text). 
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FIG. 29: (Color online) DR scattering amplitudes K r as de- 
fined in Eq. ID2I for the aa and ab processes, as a function of 
the adimensional momentum k = 2khvF / el- The real and 
imaginary part of the complex number K r are plotted as 
two different lines. Calculations are done using ez = 2.4 eV, 
7 = 84 meV, and hv F = 6.49 eVA. 



later) ££,) = 



kdk — — — , , . 

Z7T 



(D2) 



where the label a = aa or ab, and Kj^(k) are the K 
scattering amplitudes of Eqs. [DTI calculated at e L and q r , 
averaged over the angular dependence of k. 



Fig. [29] shows K^ a (k) and K^ b (k) for realistic values of 
the parameters ez, 7 and vp- Both K^ a (k) and K^ b (k) 
have a maximum near k = ez/ '(2Hvf) which corresponds 
to the DR condition (k = 1 in Fig. [29]). First we remark 
that, for realistic values of 7, the real, Re, and imaginary 
parts, Im, of the K r amplitudes arc of the same order 
of magnitude. Thus, the K r cannot be approximated 
as purely real or purely imaginary numbers. Second we 
notice that Re(-K^ h ) and lm(K^ b ) do not change their 
sign when plotted as a function of k. On the contrary, 
Re(KJJ and Im(i^) change their sign (Fig. [29]). Be- 
cause of this, the K^ b (k) inside the integral of Eq. ID2I 
add coherently, while the K^ a (k) interfere in a destruc- 
tive way. As a consequence, I u b 3> I aa , despite the fact 
that K r ab and K^ a are of the same order of magnitude. 
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